Help with syms ode solving

Hi there,

My ODE will solve when a1 = 0 but fails when a1 is not zero?

Any ideas most welcome!

function objective = myobjective(a0,a1)

syms Ca(z) Cb(z)

T=a0+a1.*z;

R= 8.3145;

Ca0 = 0.7;

k1 = 1.37.*10^10.*exp(-75000./(R.*T));

k2 = 1.19.*10^17.*exp(-125000./(R.*T));

ode1 = diff(Ca) == -4.*k1.*Ca;

ode2 = diff(Cb) == 4.*k1.*Ca - 4.*k2.*Cb;

odes = [ode1;ode2];

cond1 = Ca(0) == Ca0;

cond2 = Cb(0) == 0;

conds = [cond1; cond2];

[CaSol(z), CbSol(z)] = dsolve(odes,conds);

Cbfunction = matlabFunction(CbSol);

objective = 1./Cbfunction(2)

Thanks

Answer 1

[Star Strider]

The differential equation is nonlinear, so no analytic expression for it exists.

Try this instead —

syms Ca(z) Cb(z) z Y 
a0 = 1;
a1 = 42;
T=a0+a1.*z;
R= 8.3145;
Ca0 = 0.7;
k1 = 1.37.*10^10.*exp(-75000./(R.*T));
k2 = 1.19.*10^17.*exp(-125000./(R.*T));
ode1 = diff(Ca) == -4.*k1.*Ca;
ode2 = diff(Cb) == 4.*k1.*Ca - 4.*k2.*Cb;
odes = [ode1;ode2];
cond1 = Ca(0) == Ca0;
cond2 = Cb(0) == 0;
conds = [cond1; cond2];
% [CaSol(z), CbSol(z)] = dsolve(odes,conds)
% CaSol = vpa(CaSol,5)
% CbSol = vpa(CbSol,5)
% Cbfunction = matlabFunction(CbSol);
% objective = 1./Cbfunction(2)

[VF,Sbs] = odeToVectorField(odes);
CaCbfcn = matlabFunction(VF, 'Vars',{z,Y});

[z,CaCb] = ode15s(CaCbfcn, [0 20], [0 0.7]);

figure
plot(z,CaCb)
grid
xlabel('z')
ylabel('C')
legend(string(Sbs), 'Location','best')
set(gca, 'YLim',[min(ylim) 1])

The ‘Cb’ peak location is inversely related (with respect to ‘z’) to ‘a1’, so adjust the ‘tspan’ argument accordingly to show it correctly.

.

[Robin Thornton]

Thank you so much!

[Star Strider]

As always, my pleasure!