Coupled ODE dsolve Problem
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello! I am trying to solve the following seven ODES with the code below
syms t v(t) c(t) vxc(t) I(t) D(t) A(t) axc(t)
kf= 10;
kr= 0.15;
kf1= 3e-4;
kf2= 3e-5;
kf3= 100;
kr1= 0.5;
ode1= diff(v(t),t) == -kf*v(t)*c(t)+ kr*vxc(t) + kf2*100*I(t);
icod_1 = 0.003 ;
ode2= diff(c(t),t) == -kf*v(t)*c(t)+ kr*vxc(t) - kf3*A(t)*c(t)+ kr1*axc(t)
icod_2 = 1.0 ;
ode3= diff(vxc(t),t) == kf*v(t)*c(t)- kr*vxc(t) - kf1*I(t)
icod_3 = 0 ;
ode4= diff(I(t),t) == kf1*I(t)- kf2*I(t)
icod_4 = 0
ode5= diff(D(t),t) == kf2*I(t)
icod_5 = 0
ode6= diff(A(t),t) == - kf3*A(t)*c(t)+ kr1*axc(t)
icod_6 = 0
ode7= diff(axc(t),t) == kf3*A(t)*c(t) - kr1*axc(t)
icod_7 = 0
odes= [ode1, ode2, ode3, ode4, ode5, ode6, ode7]
IC = [icod_1,icod_2,icod_3,icod_4,icod_5,icod_6,icod_7]
[S1(t),S2(t),S3(t),S4(t),S5(t),S6(t),S7(t)] = dsolve(odes,IC)
However, I get these errors:
Error using mupadengine/feval_internal
Invalid equations.
Error in dsolve>mupadDsolve (line 334)
T = feval_internal(symengine,'symobj::dsolve',sys,x,options);
Error in dsolve (line 203)
sol = mupadDsolve(args, options);
Any suggestions to fix these errors/solve these ODE's?
0 commentaires
Réponses (1)
Star Strider
le 14 Fév 2022
Modifié(e) : Star Strider
le 15 Fév 2022
The original dsolve call output was throwing the error.
However the problem is that the system is nonlinear, and as the result the system does not have a symbolic solution (as would be the situation for most nonlinear differential equation systems).
syms t v(t) c(t) vxc(t) I(t) D(t) A(t) axc(t) Y
sympref('AbbreviateOutput',false);
kf= 10;
kr= 0.15;
kf1= 3e-4;
kf2= 3e-5;
kf3= 100;
kr1= 0.5;
ode1= diff(v(t),t) == -kf*v(t)*c(t)+ kr*vxc(t) + kf2*100*I(t)
icod_1 = v(0) == 0.003 ;
ode2= diff(c(t),t) == -kf*v(t)*c(t)+ kr*vxc(t) - kf3*A(t)*c(t)+ kr1*axc(t)
icod_2 = c(0) == 1.0 ;
ode3= diff(vxc(t),t) == kf*v(t)*c(t)- kr*vxc(t) - kf1*I(t)
icod_3 = vxc(0) == 0 ;
ode4= diff(I(t),t) == kf1*I(t)- kf2*I(t)
icod_4 = I(0) == 0;
ode5= diff(D(t),t) == kf2*I(t)
icod_5 = D(0) == 0;
ode6= diff(A(t),t) == - kf3*A(t)*c(t)+ kr1*axc(t)
icod_6 = A(0) == 0;
ode7= diff(axc(t),t) == kf3*A(t)*c(t) - kr1*axc(t)
icod_7 = axc(0) == 0;
odes= [ode1; ode2; ode3; ode4; ode5; ode6; ode7]
IC = [icod_1;icod_2;icod_3;icod_4;icod_5;icod_6;icod_7]
% [S1(t),S2(t),S3(t),S4(t),S5(t),S6(t),S7(t)] = dsolve(odes,IC)
% S = dsolve(odes,IC)
[VF,Sbs] = odeToVectorField(odes)
odes_fcn = matlabFunction(VF, 'Vars',{t,Y})
tspan = linspace(0, 1, 250);
Y0 = [1, 0.003, 0, 0, 0, 0, 0]+eps;
[t,y] = ode45(odes_fcn, tspan, Y0); % Numericaly Integrate System
figure
yyaxis left
plot(t,y(:,1))
ylabel('c')
yyaxis right
plot(t, y(:,2:end))
grid
xlabel('Time')
ylabel('All The Others')
legend(string(Sbs), 'Location','E')
The solution is to create an anonymous function of this system and integrate it numerically. The ‘odes_fcn’ can be used with ode45 or ode15s (that may be necessary if ode45 believes that this system is sufficiently stiff, and takes forever to integrate it). The ‘Sbs’ vector explains the function each equation integrates, and is useful in the legend call.
EDIT — (15 Feb 2022 at 2:15)
Added ode45 function evaluation and plot.
.
0 commentaires
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!