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

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

Answer 1

[Torsten]

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   

[Sam Chak]

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