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

조회 수: 10 (최근 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