Coupled ODE dsolve Problem

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?

Answer 1

[Star Strider]

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)

ode1 = 

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)

ode2 = 

icod_2 = c(0) == 1.0 ;
 
ode3= diff(vxc(t),t) == kf*v(t)*c(t)- kr*vxc(t) - kf1*I(t)

ode3 = 

icod_3 = vxc(0) == 0 ;
 
ode4= diff(I(t),t) == kf1*I(t)- kf2*I(t)

ode4 = 

icod_4 = I(0) == 0;
 
ode5= diff(D(t),t) == kf2*I(t)

ode5 = 

icod_5 = D(0) == 0;
 
ode6= diff(A(t),t) == - kf3*A(t)*c(t)+ kr1*axc(t)

ode6 = 

icod_6 = A(0) == 0;
 
ode7= diff(axc(t),t) ==  kf3*A(t)*c(t) - kr1*axc(t)

ode7 = 

icod_7 = axc(0) == 0;
 
odes= [ode1; ode2; ode3; ode4; ode5; ode6; ode7]

odes = 

IC = [icod_1;icod_2;icod_3;icod_4;icod_5;icod_6;icod_7]

IC = 

 
% [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)

VF = 

Sbs = 

odes_fcn = matlabFunction(VF, 'Vars',{t,Y})

odes_fcn = function_handle with value:
    @(t,Y)[Y(3).*(3.0./2.0e+1)+Y(7)./2.0-Y(1).*Y(2).*1.0e+1-Y(1).*Y(6).*1.0e+2;Y(3).*(3.0./2.0e+1)+Y(4).*(3.0./1.0e+3)-Y(1).*Y(2).*1.0e+1;Y(3).*(-3.0./2.0e+1)-Y(4).*3.0e-4+Y(1).*Y(2).*1.0e+1;Y(4).*2.7e-4;Y(4).*3.0e-5;Y(7)./2.0-Y(1).*Y(6).*1.0e+2;Y(7).*(-1.0./2.0)+Y(1).*Y(6).*1.0e+2]

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.

.