Hi all, I have never done this before, so pardon my formatting. I have a fairly long PDE solver function written, and I have attached it so as to save space. The abridged version of the problem I am having is in the pdefun subfunction of the solver. when the program reaches the first line under "%Concentration reactions," I receive the error "Attempted to access u(4); index out of bounds because numel(u)=1." There is probably a simple explanation that I am just missing. any help would be great, Thanks
function [c,f,s] = pdefun(x,t,u,DuDx)%diffusion constants:
DCO2 = 2*10^-9; %all in m^2 s^-1
DHCO3 = 1.18*10^-9;DCO3 = .995*10^-9;DH = 9.31*10^-9;DOH = 5.27*10^-9;DBOH3 = 1.11*10^-9;DBOH4 = .97*10^-9;%diffusion matrix:
Ds=[DCO2 DHCO3 DCO3 DH DOH DBOH3 DBOH4];T = 298.15; %seawater temp, degrees K
S = 35; %seawater salinity
%General Calculated Constants:
K_1 = exp(2.83655-2307.1266/T-1.5529413*log(T) ... +(-0.20760841-4.0484/T)*S^0.5+0.0846834*S-0.00654208*S^1.5); % Roy 1993
K_2 = exp(-9.226508-3351.6106/T-0.2005743*log(T) ... +(-0.106901733-23.9722/T)*S^0.5+0.1130822*S-0.00846934*S^1.5); K_w = exp(-13847.26/T+148.9652-23.6521*log(T) ... +(118.67/T-5.977+1.0495*log(T))*S^0.5-0.01615*S); % DOE 1994
K_B = exp((-8966.90-2890.53*S^.5-77.942*S+1.728*S^1.5-0.0996*S^2)/T ... +148.0248+137.1942*S^.5+1.62142*S -(24.4344+25.085*S^.5+0.2474*S)*log(T)+0.053105*(S^.5)*T); % DOE 1994%Carbonate Chemistry Reaction Constants:
kp1 = exp(1246.98-6.19e4/T-183.0*log(T)); % Schulz et al. 2006
km1 = kp1/K_1; % Schulz et al. 2006A_4 = 499002.24*exp(4.2986e-4*S^2+5.75499e-5*S); % Schulz et al. 2006kp4 = A_4*exp(-90166.83/8.31451/T)/K_w; % Schulz et al. 2006km4 = kp4*K_w/K_1; % Schulz et al. 2006kp5H = 5.0e10; % Schulz et al. 2006km5H = kp5H*K_2; % Schulz et al. 2006kp5OH = 6.0e9; % Schulz et al. 2006km5OH = kp5OH*K_w/K_2; % Schulz et al. 2006kp6 = 1.40e-3; % Schulz et al. 2006km6 = kp6/K_w; % Schulz et al. 2006km7 = 10^10;kp7 = 20;%concentration reactions:
CO2 = ((km1*u(4) + km4)*u(2) - (kp1 + kp4*u(5))*u(1))/Ds(1,1);HCO3 = (kp1*u(1) - km1*u(4)*u(2) + kp4*u(1)*u(5) - km4*u(2) + kp5H*u(4)*u(3) - km5H*u(2) + km5OH*u(3) - kp5OH*u(5)*u(2))/Ds(2,1);CO3 = (km5H*u(2) - kp5H*u(4)*u(3)+ kp5OH*u(5)*u(2) - km5OH*u(3))/Ds(3,1);H = ((km5H - km1*u(4))*u(2) + kp1*u(1) - kp5H*u(4)*u(3) + kp6 + km6*u(5)*u(4) + kp7*u(6) - km7*u(4)*u(7))/Ds(4,1);OH = (km4*u(2) - kp4*u(1)*u(5) -kp5OH*u(5)*u(2)+km5OH*u(3) + kp6 + km6*u(4)*u(5))/Ds(5,1);BOH3 = (-kp7*u(6) + km7*u(4)*u(7))/Ds(6,1);BOH4 = (kp7*u(6) - km7*u(4)*u(7))/Ds(7,1);%Set PDE Values
c=ones(1,7)/Ds;f=ones(1,7) .*DuDx;s=[CO2; HCO3; CO3; H; OH; BOH3; BOH4];
Best Answer