Solving system of odes with a power using ode45

I have the following system of first order ode i would like to solve it using ode45 1) dX/dt = -0.000038*X - (X*(X/Xinit)^frac)*rext 2) dY/dt = - 0.000038*Y + rext*X - rtra*Y + Sr 3) dZ/dt = - 0.000038*Z + rext*Y - rtra*Z + Sti 4) dU/dt = 0.000038*U + rext*Z - rvol*U + Sfeu Satisfying X(0)=Y(0)=Z(0)=U(0)=0 Where the functions are X, Y,Z and U and the variable is t. The others parameters are known constant It is possible toi solve it with ode45 ? Since a power appear in the first equation

Answer 1

[Sam Chak]

Hi @Thomas TJOCK-MBAGA

I presume that Xinit is not equal to , and . I tested this with the ode45 solver, and it also works with non-zero initial values.

Method 1: Using the function ode45() directly

% Parameters
Xinit = 1;
frac  = 1/2;
rext  = 1;
rtra  = 1;
rvol  = 1;
Sr    = 1;
Sti   = 1;
Sfeu  = 1;

% Create a function handle F for a system of 1st-order ODEs:
F      = @(t, y) [- 0.000038*y(1) - (y(1)*(y(1)/Xinit)^frac)*rext;
                  - 0.000038*y(2) + rext*y(1) - rtra*y(2) + Sr;
                  - 0.000038*y(3) + rext*y(2) - rtra*y(3) + Sti;
                    0.000038*y(4) + rext*y(3) - rvol*y(4) + Sfeu];
tspan  = [0 10];
y0     = [3; 2; 1; 0];
[t, y] = ode45(F, tspan, y0);

% Plotting the result
plot(t, y, "-o"), grid on
xlabel('Time, t (seconds)'), ylabel('\bf{y}(t)')
legend('X', 'Y', 'Z', 'U', 'location', 'NW')

Method 2: Using the 'ode' object (introduced in R2023b)

% Parameters
Xinit = 1;
frac  = 1/2;
rext  = 1;
rtra  = 1;
rvol  = 1;
Sr    = 1;
Sti   = 1;
Sfeu  = 1;

% Setting up the ODE object:
F              = ode;
F.InitialValue = [3; 2; 1; 0];
F.ODEFcn       = @(t, y) [- 0.000038*y(1) - (y(1)*(y(1)/Xinit)^frac)*rext;
                          - 0.000038*y(2) + rext*y(1) - rtra*y(2) + Sr;
                          - 0.000038*y(3) + rext*y(2) - rtra*y(3) + Sti;
                            0.000038*y(4) + rext*y(3) - rvol*y(4) + Sfeu];
F.Solver       = "ode45";
S              = solve(F, 0, 10);   % Solving F over the time from 0 to 10 s

% Plotting the result
plot(S.Time, S.Solution, "-o"), grid on
xlabel('Time, t (seconds)'), ylabel('\bf{y}(t)')
legend('X', 'Y', 'Z', 'U', 'location', 'NW')

[William Rose]

@Sam Chak, nice!

[Sam Chak]

@William Rose, Thanks. I'm still testing the newly introduced ode object in R2023b. Unfortunately, the ode object solver method does not permit the implicit ODE solver ode15i() at the moment.

[William Rose]

@Thomas TJOCK-MBAGA ,

You said "...Satisfying X(0)=Y(0)=Z(0)=U(0)=0".

You may want to adjust the solution of @Sam Chak, to use

F.InitialValue = [0; 0; 0; 0];

instead of

F.InitialValue = [3; 2; 1; 0];

@Sam Chak, I have not learned about using "solver" the way you have done to set up and solve ODEs. That is good to know about. Thanks.

[Thomas TJOCK-MBAGA]

Thanks you i try it. This program also work in MATLAB 2016a?

[Thomas TJOCK-MBAGA]

I try the code with MATLAB 2016a an error occurs indefined variable or function "ode"

[Steven Lord]

No, this code will not work in release R2016a. The ode object this code uses was introduced in release R2023b.

[Thomas TJOCK-MBAGA]

Which modification i am suppose to do??

[William Rose]

@Thomas TJOCK-MBAGA,

To run the code in R2016, use the code ideas of @Sam Chak, translated to the format described in the examples which you will find in the help for ode45().

[Sam Chak]

Hi @Thomas TJOCK-MBAGA

I have updated my Answer so that you can use Method #1 ODE solver in your MATLAB R2016a version. If you find the solution helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Thanks a bunch! By the way, you should also vote for other helpful answers as tokens of appreciation.

[Thomas TJOCK-MBAGA]

Thanks you, it work very well

[Thomas TJOCK-MBAGA]

Non i have the same system but with thd first equation be a PDE i would like to solve it pdepe solver and using ode45 1) 1) R*dX/dt = -0.000038*X - (X*(X/Xinit)^frac)*rext -v*dX/dx + alpha*d^2X/dx^2 2) dY/dt = - 0.000038*Y + rext*X - rtra*Y + Sr 3) dZ/dt = - 0.000038*Z + rext*Y - rtra*Z + Sti 4) dU/dt = 0.000038*U + rext*Z - rvol*U + Sfeu Satisfying Y(0)=Z(0)=U(0)=0 X(t,0)=0; dX/dx(t,x=L)=0; X(t=0,x)=Xinit Where the functions are X, Y,Z and U and the variable are x and t. The others parameters are known constant. If frac=0 i know how to.solve it with Laplace transform and then by an intégration over the space domaine i obtain X(t) and then the functions. No

[William Rose]

@Thomas TJOCK-MBAGA, sounds like a good new question.

Answer 2

[William Rose]

@Thomas TJOCK-MBAGA,

Yes you can do it with ode45().

dX/dt = -0.000038*X - (X*(X/Xinit)^frac)*rext

dY/dt = - 0.000038*Y + rext*X - rtra*Y + Sr

dZ/dt = - 0.000038*Z + rext*Y - rtra*Z + Sti

dU/dt = 0.000038*U + rext*Z - rvol*U + Sfeu

You want to replace X,Y,Z,U with x(1),x(2),x(3),x(4).

Your equation for dX/dt includes (X/Xinit)^frac. If Xinit=X(0), you have a divide-by-zero problem, since you said X(0)=0.

[Thomas TJOCK-MBAGA]

Thanks

[Thomas TJOCK-MBAGA]

Xinit IS different to 0 Thanks