system of linear differential equation with Runge Kutta Method 4 order

dT/dt= aT(1-bT) - (cN +jD+kL)T - KzT

[Diksha Gautam]

Hello sir,

please check my code and resolve the error

mu = 0.012;

lambda = 1.2e-6;

omega = 0.75;

rho = 0.1;

K_T = 4e-4;

upsilon = 2.5e-4;

M_mt = 0.1;

alpha = 2.5e-3;

K_N = 4e-5;

epsilon = 1.2e-3;

beta = 0.3;

eta = 0.02;

sigma = 0.02;

K_D = 2e-4;

j = 4e-4;

gamma = 0.06;

P_LI = 2e-4;

K_L = 4e-4;

g = 1.2e-4;

v_l = 0.01;

% Time horizon

t0 = 0; tf = 12;

% Initial conditions

T0 = 1; N0 = 1e-5; D0 = 1e-5; L0 = 1e-5;

% Time delays

tau4 = 1; tau5 = 1; tau6 = 1; tau7 = 1;

% Solve the optimal control problem

x0 = [T0, N0, D0, L0];

[t, x] = ode45(@(t,x) ode_cancer(t, x, tau4, tau5, tau6, tau7, lambda, omega, rho, K_T, upsilon, M_mt, alpha, K_N, epsilon, beta, eta, sigma, K_D, j, gamma, P_LI, K_L, g, v_l), [t0, tf], x0);

% Extract the states

T = x(:,1);

N = x(:,2);

D = x(:,3);

L = x(:,4);

% Compute the optimal therapy function

Lambda = zeros(size(t));

for i = 1:length(t)

pT = -lambda*(1-mu*T(i)) + (omega*N(i) + rho*D(i) + gamma*L(i)) - K_T*Lambda(i)*T(i);

pN = upsilon + M_mt - (omega*T(i) - alpha*D(i))*N(i) - K_N*Lambda(i)*N(i) - epsilon*N(i);

pD = beta - (eta*L(i) + alpha*N(i) - sigma*T(i))*D(i) - K_D*Lambda(i)*D(i) - j*D(i);

pL = rho*D(i)*T(max(1,i-tau2))*T(max(1,i-tau4)) - gamma*L(i) - chi*N(i)*L(max(1,i-tau3))^2 + omega*N(i)*T(max(1,i-tau1))*T(max(1,i-tau4)) + P_LI - K_L*Lambda(i)*L(i) - g*L(i);

H = pT*Lambda(i) + pN*Lambda(i) + pD*Lambda(i) + pL*Lambda(i);

if H > 0

Lambda(i) = 1;

else

Lambda(i) = 0;

end

end

% Plot the results

figure;

subplot(2,2,1);

plot(t, T);

xlabel('Time (days)');

ylabel('Tumor Size (cm^3)');

title('Tumor Growth');

subplot(2,2,2);

plot(t, Lambda);

xlabel('Time (days)');

ylabel('Treatment');

title('Optimal Treatment');

[Torsten]

Hello sir,

please check my code and resolve the error

Not possible without your ode function.

[Sam Chak]

Hi @Diksha Gautam,

Please show the code of the ode_cancer() function.

[t, x] = ode45(@(t,x) ode_cancer(t, x, tau4, tau5, tau6, tau7, lambda, omega, rho, K_T, upsilon, M_mt, alpha, K_N, epsilon, beta, eta, sigma, K_D, j, gamma, P_LI, K_L, g, v_l), [t0, tf], x0);

Unrecognized function or variable 'ode_cancer'.

Error in solution>@(t,x)ode_cancer(t,x,tau4,tau5,tau6,tau7,lambda,omega,rho,K_T,upsilon,M_mt,alpha,K_N,epsilon,beta,eta,sigma,K_D,j,gamma,P_LI,K_L,g,v_l) (line 29)
[t, x] = ode45(@(t,x) ode_cancer(t, x, tau4, tau5, tau6, tau7, lambda, omega, rho, K_T, upsilon, M_mt, alpha, K_N, epsilon, beta, eta, sigma, K_D, j, gamma, P_LI, K_L, g, v_l), [t0, tf], x0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Answer 1

[Sam Chak]

@Diksha Gautam

You can follow the examples to use ode45():

https://www.mathworks.com/help/matlab/ref/ode45.html

% constants
s = [1.3e4 3.3e4 5.3e4 7.3e4];
C = 7.13e-10;
p = 5e-2;
D = 1;
q = 4.23e3;
L = 1;
g = 2.4e-2;
h = 1e4;
T = 1;
z = 6e-1;
M = 1e4;
e = 4.12e-2;

for j = 1:length(s)
    % ODE
    f = @(t, N) s(j) + (g*N*T^2)/(h + T^2) - (C*T - p*D - q*L)*N - z*M*N - e*N;
    tspan = [0 0.01];
    init = 200;

    % ODE45 solver
    [t, N] = ode45(f, tspan, init);
    plot(t, N, 'linewidth', 1.5)
    hold on
end

hold off
grid on
xlabel('t')
ylabel('N(t)')

[Diksha Gautam]

Thank you.. But I applied this method to other equation then got error

[Sam Chak]

No issue on mine.

clear all; clc

% constants
s = 1.3e4;
C = 7.13e-10;
p = 5e-2;
D = 1;
q = 4.23e3;
L = 1;
g = 2.4e-2;
h = 1e4;
T = 1;
z = 6e-1;
M = 1e4;
e = 4.12e-2;

% Ordinary differential equations
f = @(t, N) s + (g*N*T^2)/(h + T^2) - (C*T - p*D - q*L)*N - z*M*N - e*N;
tspan = [0 0.01];
init = 200;

% Runge-Kutta Dormand-Prince 4/5 solver
[t, N] = ode45(f, tspan, init);

plot(t, N, 'linewidth', 1.5)

[Diksha Gautam]

Thank you so much

[Diksha Gautam]

If I have different values of s then want to plot all values in one graph by different colors so how can i will do?

Answer 2

[Diksha Gautam]

Can you please tell me how can i plot multiple graph with function and if possible take above example

[Torsten]

It's unclear what your problem is.

If f and g are functions and x is a column vector,

plot (x,[f(x),g(x)])

plots f(x) and g(x) against x.

[Diksha Gautam]

hello sir please provide code for hamiltonian equation to run and plot graph with default parameters

[Diksha Gautam]

@Sam Chak please share the code for

ddtΤt= λΤt1-μΤt-ωNt-τ1+ρDt-τ2+γLt-τ3Τt-ΚΤΛtΤt-τ4

[Sam Chak]

Sorry, I'm outside. Didn't the original problem solved 10 months ago. What exactly is this? I cannot interpret it mathematically in my mind.

Answer 3

[Sam Chak]

Hi @Diksha Gautam

I have updated the code for running a few iterations to generate the results for different values of s.

Please consider voting 👍 this for support.

[Diksha Gautam]

Ok..thank you so much

[Diksha Gautam]

Hii @Sam Chak Can you please tell me for rkm5 order which function is apply in matlab Is ode54?

[Sam Chak]

Hi @Diksha Gautam,

ode54() does not exist. But ode45() is based on an explicit Runge-Kutta (4th- and 5th-order) formula, known as the Dormand-Prince method.

help ode54

--- ode54 not found. Showing help for ode45 instead. ---

 ODE45  Solve non-stiff differential equations, medium order method.
    [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0) integrates the system of
    differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end)
    with initial conditions Y0. Each row in the solution array YOUT
    corresponds to a time in the column vector TOUT. 
      * ODEFUN is a function handle. For a scalar T and a vector Y,
        ODEFUN(T,Y) must return a column vector corresponding to f(t,y).
      * TSPAN is a two-element vector [T0 TFINAL] or a vector with
        several time points [T0 T1 ... TFINAL]. If you specify more than
        two time points, ODE45 returns interpolated solutions at the
        requested times.
      * YO is a column vector of initial conditions, one for each equation.
 
    [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration
    option values in the fields of a structure, OPTIONS. Create the
    options structure with odeset.
 
    [TOUT,YOUT,TE,YE,IE] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) produces
    additional outputs for events. An event occurs when a specified function
    of T and Y is equal to zero. See ODE Event Location for details.
 
    SOL = ODE45(...) returns a solution structure instead of numeric
    vectors. Use SOL as an input to DEVAL to evaluate the solution at
    specific points. Use it as an input to ODEXTEND to extend the
    integration interval.
 
    ODE45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
    nonsingular. Use ODESET to set the 'Mass' property to a function handle
    or the value of the mass matrix. ODE15S and ODE23T can solve problems
    with singular mass matrices.
 
    ODE23, ODE45, ODE78, and ODE89 are all single-step solvers that use
    explicit Runge-Kutta formulas of different orders to estimate the error
    in each step.
      * ODE45 is for general use.
      * ODE23 is useful for moderately stiff problems.
      * ODE78 and ODE89 may be more efficient than ODE45 on non-stiff problems
        that are smooth except possibly for a few isolated discontinuities.
      * ODE89 may be more efficient than ODE78 on very smooth problems, when 
        integrating over long time intervals, or when tolerances are tight.
 
    Example
          [t,y]=ode45(@vdp1,[0 20],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution. 
    
    Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
      float: double, single
 
    See also ODE23, ODE78, ODE89, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB,
             ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, FUNCTION_HANDLE.

    Documentation for ode45
       doc ode45

    Other uses of ode45

       dlarray/ode45

[Diksha Gautam]

Ok. Thank you but In matlab any function to solve higher order equation by rkm such 5 n 6 order?

[Sam Chak]

Hi @Diksha Gautam

If you looking for the "most efficient", try ode78()

help ode78

 ODE78  Solve non-stiff differential equations, high order method.
    [TOUT,YOUT] = ODE78(ODEFUN,TSPAN,Y0) integrates the system of
    differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end)
    with initial conditions Y0. Each row in the solution array YOUT
    corresponds to a time in the column vector TOUT.
      * ODEFUN is a function handle. For a scalar T and a vector Y,
        ODEFUN(T,Y) must return a column vector corresponding to f(t,y).
      * TSPAN is a two-element vector [T0 TFINAL] or a vector with
        several time points [T0 T1 ... TFINAL]. If you specify more than
        two time points, ODE78 returns interpolated solutions at the
        requested times.
      * YO is a column vector of initial conditions, one for each equation.
 
    [TOUT,YOUT] = ODE78(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration
    option values in the fields of a structure, OPTIONS. Create the
    options structure with odeset.
 
    [TOUT,YOUT,TE,YE,IE] = ODE78(ODEFUN,TSPAN,Y0,OPTIONS) produces
    additional outputs for events. An event occurs when a specified function
    of T and Y is equal to zero. See ODE Event Location for details.
 
    SOL = ODE78(...) returns a solution structure instead of numeric
    vectors. Use SOL as an input to DEVAL to evaluate the solution at
    specific points. Use it as an input to ODEXTEND to extend the
    integration interval.
 
    ODE78 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
    nonsingular. Use ODESET to set the 'Mass' property to a function handle
    or the value of the mass matrix. ODE15S and ODE23T can solve problems
    with singular mass matrices.
 
    ODE23, ODE45, ODE78, and ODE89 are all single-step solvers that use
    explicit Runge-Kutta formulas of different orders to estimate the error
    in each step.
      * ODE45 is for general use.
      * ODE23 is useful for moderately stiff problems.
      * ODE78 and ODE89 may be more efficient than ODE45 on non-stiff problems
        that are smooth except possibly for a few isolated discontinuities.
      * ODE89 may be more efficient than ODE78 on very smooth problems, when
        integrating over long time intervals, or when tolerances are tight.
      * ODE78 and ODE89 may not be as fast or as accurate as ODE45 in single
        precision.
 
    Example
          opts = odeset('AbsTol',1e-8,'RelTol',1e-6);
          [t,y] = ode78(@vdp1,[0 20],[2 0],opts);
          plot(t,y(:,1));
      solves the system y' = vdp1(t,y), using the specified tolerances, and
      plots the first component of the solution.
 
    Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
      float: double, single
 
    See also ODE89, ODE45, ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB,
             ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, FUNCTION_HANDLE.

    Documentation for ode78
       doc ode78

If you looking for the "most robust", try ode89()

help ode89

 ODE89  Solve non-stiff differential equations, high order method.
    [TOUT,YOUT] = ODE89(ODEFUN,TSPAN,Y0) integrates the system of
    differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end)
    with initial conditions Y0. Each row in the solution array YOUT
    corresponds to a time in the column vector TOUT.
      * ODEFUN is a function handle. For a scalar T and a vector Y,
        ODEFUN(T,Y) must return a column vector corresponding to f(t,y).
      * TSPAN is a two-element vector [T0 TFINAL] or a vector with
        several time points [T0 T1 ... TFINAL]. If you specify more than
        two time points, ODE89 returns interpolated solutions at the
        requested times.
      * YO is a column vector of initial conditions, one for each equation.
 
    [TOUT,YOUT] = ODE89(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration
    option values in the fields of a structure, OPTIONS. Create the
    options structure with odeset.
 
    [TOUT,YOUT,TE,YE,IE] = ODE89(ODEFUN,TSPAN,Y0,OPTIONS) produces
    additional outputs for events. An event occurs when a specified function
    of T and Y is equal to zero. See ODE Event Location for details.
 
    SOL = ODE89(...) returns a solution structure instead of numeric
    vectors. Use SOL as an input to DEVAL to evaluate the solution at
    specific points. Use it as an input to ODEXTEND to extend the
    integration interval.
 
    ODE89 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
    nonsingular. Use ODESET to set the 'Mass' property to a function handle
    or the value of the mass matrix. ODE15S and ODE23T can solve problems
    with singular mass matrices.
 
    ODE23, ODE45, ODE78, and ODE89 are all single-step solvers that use
    explicit Runge-Kutta formulas of different orders to estimate the error
    in each step.
      * ODE45 is for general use.
      * ODE23 is useful for moderately stiff problems.
      * ODE78 and ODE89 may be more efficient than ODE45 on non-stiff problems
        that are smooth except possibly for a few isolated discontinuities.
      * ODE89 may be more efficient than ODE78 on very smooth problems, when
        integrating over long time intervals, or when tolerances are tight.
      * ODE78 and ODE89 may not be as fast or as accurate as ODE45 in single
        precision.
 
    Example
          opts = odeset('AbsTol',1e-10,'RelTol',1e-8);
          [t,y] = ode89(@vdp1,[0 20],[2 0],opts);
          plot(t,y(:,1));
      solves the system y' = vdp1(t,y), using the specified tolerances, and
      plots the first component of the solution.
 
    Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
      float: double, single
 
    See also ODE78, ODE45, ODE23, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB,
             ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, FUNCTION_HANDLE.

    Documentation for ode89
       doc ode89

Answer 4

[Diksha Gautam]

Thank you

Answer 5

[Diksha Gautam]

function [T, N, D, L, PT, PN, PD, PL, Lambda] = ode_optimization()

% Define the system parameters

mu = 0.1; omega = 0.2; rho = 0.15; gamma = 0.1; K_T = 0.2; K_N = 0.2;

K_D = 0.2; K_L = 0.2; eta = 0.1; alpha = 0.1; sigma = 0.1; epsilon = 0.1;

chi = 0.1; delta = 0.1; v_l = 0.1; P_LI = 0.1; M_mt = 0.1; beta = 0.1;

tau_1 = 0.2; tau_2 = 0.3; tau_3 = 0.4; tau_4 = 0.5; tau_5 = 0.6; tau_6 = 0.7;

tau_7 = 0.8; lambda_min = 0.0012; lambda_max = 1;

% Define the time interval and initial conditions

tspan = [0 10];

T0 = 0.1; N0 = 0.1; D0 = 0.1; L0 = 0.1;

PT0 = 0; PN0 = 0; PD0 = 0; PL0 = 0;

Lambda0 = 0;

% Define the objective function

J = @(T) T(end);

% Define the Hamiltonian

H = @(t, T, N, D, L, PT, PN, PD, PL, Lambda, lambda) ...

PT*(lambda*T*(1-mu*T) - (omega*N + rho*D + gamma*L)*T - K_T*Lambda*T(t-tau_4)) ...

+ PN*(v_l + M_mt - (omega*T - alpha*D)*N - K_N*Lambda*N(t-tau_5) - epsilon*N) ...

+ PD*(beta - (eta*L + alpha*N - sigma*T)*D - K_D*Lambda*D(t-tau_6) - delta*D) ...

+ PL*(rho*D*T(t-tau_2)*T(t-tau_4) - gamma*L - chi*N(t-tau_5)*L(t-tau_3)^2 + omega*N(t-tau_1)*T(t-tau_4) + P_LI ...

- K_L*Lambda*L(t-tau_7) - delta*L + v_l);

% Define the ODE system

ode_system = @(t, y, lambda) [H(t, y(1), y(2), y(3), y(4), y(5), y(6), y(7), y(8), y(9), lambda); % H(T)

y(1)*(1-y(1)) - eta*y(4)*y(2); % dT/dt

alpha*y(1)*y(3) - omega*y(2); % dN/dt

sigma*y(1)*y(3) - eta*y(4)*y(2); % dD/dt

chi*y(5)*y(3)^2 - gamma*y(4) - delta*y(4) + omega*y(2)*y(2-tau_1)*y(1-tau_4) + P_LI - K_L*y(5-tau_7)*Lambda - v_l*y(4); % dL/dt

-PT*(lambda*T*(1-mu*T) - (omega*N + rho*D + gamma*L)*T - K_T*Lambda*T(t-tau_4)); % dPT/dt

-PN*(v_l + M_mt - (omega*T - alpha*D)*N - K_N*Lambda*N(t-tau_5) - epsilon*N); % dPN/dt

-PD*(beta - (eta*L + alpha*N - sigma*T)*D - K_D*Lambda*D(t-tau_6) - delta*D); % dPD/dt

-PL*(rho*D*T(t-tau_2)*T(t-tau_4) - gamma*L - chi*N(t-tau_5)*L(t-tau_3)^2 + omega*N(t-tau_1)*T(t-tau_4) + P_LI - K_L*Lambda*L(t-tau_7) - delta*L + v_l)]; % dLambda/dt

% Solve the ODE system

sol = dde23(@(t,y,Z,lambda) ode_system(t,y,lambda), ...

[tau_7 tau_6 tau_5 tau_4 tau_3 tau_2 tau_1], ...

[L0; D0; N0; T0; PL0; PD0; PN0; PT0; Lambda0], ...

tspan, ...

@(t) interp1(sol.x, sol.y, t));

% Plot the solution

figure

plot(sol.x, sol.y(1,:), 'r', ...

sol.x, sol.y(2,:), 'g', ...

sol.x, sol.y(3,:), 'b', ...

sol.x, sol.y(4,:), 'm')

legend('L', 'D', 'N', 'T')

xlabel('Time')

ylabel('optimization')

please check code has error how i will resolve it

Answer 6

[Diksha Gautam]

hii @Sam Chak please check this code

% Define the right-hand side of the DDE

dde = @(t,y,Z) dde_rhs(t,y,Z,lambda_val,mu_val,omega_val,rho_val,gamma_val,K_T_val,Lambda_val,tau_1,tau_2,tau_3,tau_4);

% Set the parameter values and delays

lambda_val = 0.5;

mu_val = 0.1;

omega_val = 0.2;

rho_val = 0.3;

gamma_val = 0.4;

K_T_val = 0.6;

Lambda_val = 0.7;

tau_1 = 1.0;

tau_2 = 2.0;

tau_3 = 3.0;

tau_4 = 4.0;

T_end = 10.0;

% Set the initial conditions

T0 = 1.0;

T_init = @(t) T0*ones(size(t));

Z_init = [T_init(tau_4:-tau_1:0), zeros(3,length(tau_1:0))];

% Set the time grid and options for dde23 solver

tspan = [tau_4, T_end];

options = ddeset('RelTol',1e-6,'AbsTol',1e-6);

% Solve the DDE using dde23 solver

sol = dde23(dde,[tau_1,tau_2,tau_3,tau_4],Z_init,tspan,options);

% Get the solution for T(t)

T_sol = sol.y(1,:);

% Plot the solution

figure;

plot(sol.x, T_sol, 'LineWidth', 2);

xlabel('Time');

ylabel('T(t)');

title('Solution of Delay Differential Equation');

grid on