Matlab numerical simulation for ode system by changing the parameter value in range graph not running error coming
Show older comments
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);
Accepted Answer
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()
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 Comment
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.
), 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,:)')
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
More Answers (0)
Categories
Find more on Guidance, Navigation, and Control (GNC) in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
