Matlab numerical simulation for ode system by changing the parameter value in range graph not running error coming

5 ビュー (過去 30 日間)
Dhivyadharshini
Dhivyadharshini 2025 年 9 月 3 日
コメント済み: Sam Chak 2025 年 9 月 3 日
function su
options = odeset('RelTol',1e-6,'Stats','on');
%initial conditions
Xo = [0.005; 0.0007; 0.001; 0.0001; 0.001; 0; 0; 0.007];
tspan = [0,120];
tic
[t,X] = ode45(@TestFunction,tspan,Xo,options);
toc
%figure
%plot(t, X(:,1), 'b')
plot(t, X(:,2), 'b')
%plot(t, X(:,3), 'g')
%plot(t, X(:,4), 'm')
%plot(t, X(:,5), 'k')
%plot(t, X(:,6), 'c')
%plot(t, X(:,7), 'y')
%plot(t, X(:,8), 'y')
hold on
% legend('x1','x2')
% ylabel('x - Population')
% xlabel('t - Time')
hold on
return
function [dx_dt]= TestFunction(~,x)
% Parameters
s = 0.038;
alpha = 0.02;
gamma = 0.10;
r = 0.03;
dH = 0.0083;
dX = 0.0125;
rho = 0.07;
K = 500;
beta = 0.0005;
theta = 0.03;
eta = 0.015;
muH = 0.015;
muI = 0.08;
UP = 0.20;
deltaP = 0.033;
UL = 0.50;
deltaL = 0.05;
wP = 10000;
wL = 20000;
wV = 1;
LC = 0.05;
dx_dt(1) = s - (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - dH * x(1) + r * x(2);
dx_dt(2) = (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - (r+dX) * x(2);
dx_dt(3) = rho * x(3) * (1 - (x(3) + x(4)) / K) - beta * x(3) * x(5) - muH * x(3) - LC * x(8) * x(3);
dx_dt(4) = beta * x(3) * x(5) - muI * x(4) - LC * x(8) * x(4);
dx_dt(5) = theta * x(4) - eta * x(5) - beta * x(3) * x(5);
dx_dt(6) = UP - deltaP * x(6);
dx_dt(7) = UL - deltaL * x(7);
dx_dt(8) = wP * x(6) + wL * x(7) + wV * x(5);
% dx_dt(1) = (s) - ((alpha1 * (x(3)+x(4)) * x(1)) / (1 + gamma1 * x(1))) - ((alpha2 * x(5) * x(1)) / (1 + gamma2 * x(1)));
%dx_dt(2) = ((alpha1 * (x(3)+x(4)) * x(1)) / (1 + gamma1 * x(1)) + (alpha2 * x(5) * x(1)) / (1 + gamma2 * x(1))) - ((k1 * x(3) + k2 * x(5)) * x(2)) - (r * x(2));
%dx_dt(3) = (gamma1 * x(3) * (1 - ((x(3)+x(4)) / kR))) - (muR * x(3)) - (beta * x(3) * x(4));
%dx_dt(4) = (beta * x(3) * x(4)) - ((muR + gamma * x(6)) * x(4));
%dx_dt(5) = (gamma2 * x(5) * (1 - x(5) / kW)) - (muW * x(5)) - (theta * (x(3)+x(4)) * x(5));
% dx_dt(6) = (delta_v) + (rho * x(4)) - (eta * x(6));
dx_dt = dx_dt';
return
hold on In this code error coming. This is the error.
Operation terminated by user during ode45
In summ (line 7)
[t,X] = ode45(@TestFunction,tspan,Xo,options);

サインインしてこの質問に回答する。

採用された回答

Torsten
Torsten 2025 年 9 月 3 日
編集済み: Torsten 2025 年 9 月 3 日
Ran in:
Your ODE system is stiff. Use "ode15s" instead of "ode45".
The reason for the error you receive seems to be that you interrupted the computation (maybe because it took too long).
su()
121 successful steps 29 failed attempts 575 function evaluations 27 partial derivatives 52 LU decompositions 310 solutions of linear systems Elapsed time is 0.151195 seconds.
function su
options = odeset('RelTol',1e-6,'Stats','on');
%initial conditions
Xo = [0.005; 0.0007; 0.001; 0.0001; 0.001; 0; 0; 0.007];
tspan = [0,120];
tic
[t,X] = ode15s(@TestFunction,tspan,Xo,options);
toc
%figure
%plot(t, X(:,1), 'b')
plot(t, X(:,2), 'b')
%plot(t, X(:,3), 'g')
%plot(t, X(:,4), 'm')
%plot(t, X(:,5), 'k')
%plot(t, X(:,6), 'c')
%plot(t, X(:,7), 'y')
%plot(t, X(:,8), 'y')
%hold on
% legend('x1','x2')
% ylabel('x - Population')
% xlabel('t - Time')
%hold on
end
function [dx_dt]= TestFunction(~,x)
% Parameters
s = 0.038;
alpha = 0.02;
gamma = 0.10;
r = 0.03;
dH = 0.0083;
dX = 0.0125;
rho = 0.07;
K = 500;
beta = 0.0005;
theta = 0.03;
eta = 0.015;
muH = 0.015;
muI = 0.08;
UP = 0.20;
deltaP = 0.033;
UL = 0.50;
deltaL = 0.05;
wP = 10000;
wL = 20000;
wV = 1;
LC = 0.05;
dx_dt(1) = s - (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - dH * x(1) + r * x(2);
dx_dt(2) = (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - (r+dX) * x(2);
dx_dt(3) = rho * x(3) * (1 - (x(3) + x(4)) / K) - beta * x(3) * x(5) - muH * x(3) - LC * x(8) * x(3);
dx_dt(4) = beta * x(3) * x(5) - muI * x(4) - LC * x(8) * x(4);
dx_dt(5) = theta * x(4) - eta * x(5) - beta * x(3) * x(5);
dx_dt(6) = UP - deltaP * x(6);
dx_dt(7) = UL - deltaL * x(7);
dx_dt(8) = wP * x(6) + wL * x(7) + wV * x(5);
% dx_dt(1) = (s) - ((alpha1 * (x(3)+x(4)) * x(1)) / (1 + gamma1 * x(1))) - ((alpha2 * x(5) * x(1)) / (1 + gamma2 * x(1)));
%dx_dt(2) = ((alpha1 * (x(3)+x(4)) * x(1)) / (1 + gamma1 * x(1)) + (alpha2 * x(5) * x(1)) / (1 + gamma2 * x(1))) - ((k1 * x(3) + k2 * x(5)) * x(2)) - (r * x(2));
%dx_dt(3) = (gamma1 * x(3) * (1 - ((x(3)+x(4)) / kR))) - (muR * x(3)) - (beta * x(3) * x(4));
%dx_dt(4) = (beta * x(3) * x(4)) - ((muR + gamma * x(6)) * x(4));
%dx_dt(5) = (gamma2 * x(5) * (1 - x(5) / kW)) - (muW * x(5)) - (theta * (x(3)+x(4)) * x(5));
% dx_dt(6) = (delta_v) + (rho * x(4)) - (eta * x(6));
dx_dt = dx_dt';
end
  1 件のコメント
Sam Chak
Sam Chak 2025 年 9 月 3 日
Ran in:
Hi @Dhivyadharshini
For your records, state no.8 grows rapidly to become very large (at the magnitude of ), which may contribute to the system becoming increasingly stiff over time.
options = odeset('RelTol', 1e-6);
Xo = [0.005; 0.0007; 0.001; 0.0001; 0.001; 0; 0; 0.007];
tspan = [0, 120];
[t, X] = ode15s(@TestFunction, tspan, Xo, options);
figure
plot(t, X(:,8), 'b')
ylabel('x - Population')
xlabel('t - Time')
dx_dt = TestFunction(t(end), X(end,:)')
dx_dt = 8×1
1.0e+05 * 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 2.5895
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function [dx_dt]= TestFunction(t, x)
% Parameters
s = 0.038;
alpha = 0.02;
gamma = 0.10;
r = 0.03;
dH = 0.0083;
dX = 0.0125;
rho = 0.07;
K = 500;
beta = 0.0005;
theta = 0.03;
eta = 0.015;
muH = 0.015;
muI = 0.08;
UP = 0.20;
deltaP = 0.033;
UL = 0.50;
deltaL = 0.05;
wP = 10000;
wL = 20000;
wV = 1;
LC = 0.05;
dx_dt(1)= s - (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - dH * x(1) + r * x(2);
dx_dt(2)= (alpha * x(3) * x(1)) / (1 + gamma * x(1)) - (r+dX) * x(2);
dx_dt(3)= rho * x(3) * (1 - (x(3) + x(4)) / K) - beta * x(3) * x(5) - muH * x(3) - LC * x(8) * x(3);
dx_dt(4)= beta * x(3) * x(5) - muI * x(4) - LC * x(8) * x(4);
dx_dt(5)= theta * x(4) - eta * x(5) - beta * x(3) * x(5);
dx_dt(6)= UP - deltaP * x(6);
dx_dt(7)= UL - deltaL * x(7);
dx_dt(8)= wP * x(6) + wL * x(7) + wV * x(5);
dx_dt = dx_dt';
end

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File Exchange で Guidance, Navigation, and Control (GNC) についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by