syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H
eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;
eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;
eqn3 = FF == (d(x) * b(x)) / a(x);
eqn4 = G == (d(x) * c(x)) / b(x);
eqn5 = H == d(x) * e(x);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];
vars = [a(x); b(x); c(x); d(x); e(x)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H);
A = 2.89e-9;
B = 2.35e-9;
C = 1.69e-9;
D = 1.6e-8;
E = 9.25e-9;
FF = 6.24;
G = 5.68e-5;
H = 5.3e-8;
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H);
vars;
y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10^(-3), 'AbsTol' , 10^(-3));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')
Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]
Error in WaterProfiles>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H)
Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});
Error in WaterProfiles (line 45)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
