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
采纳的回答
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.
.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Numeric Solvers 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
