Solving system of odes with a power using ode45

3 次查看（过去 30 天）

Thomas TJOCK-MBAGA 2023-9-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2023162-solving-system-of-odes-with-a-power-using-ode45

评论： William Rose 2023-9-19

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sam Chak 2023-9-19

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2023162-solving-system-of-odes-with-a-power-using-ode45#answer_1313452

编辑：Sam Chak 2023-9-19

在 MATLAB Online 中打开

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')

12 个评论
显示 10更早的评论隐藏 10更早的评论

Thomas TJOCK-MBAGA 2023-9-19

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