How to solve ODEs with time dependent parameters by RK4 method?

Hello, everyone! I have an issue with solving ODEs by the RK4 method. The bacteria growth can be described by the couples of equation, the total number N can be divided into two parts (N1 and N2, which represented by (1) and (2)). The growth rate umax and m depend on temperature (T). T changed with time t linearly between two points (time, temp). I know there is something wrong with my code, but I have just begun to learn matlab. Does anybody can help me figure it out? Any help will be gratefully appreciated.

t = [0, 5, 10, 15, 20, 25, 30, 35, 50, 75, 100]; %time

T = [10, 4, 2, 15, 15, 2, 10, 2, 2, 15, 10]; %Temp

function main
clear all
clc
close all
t_in = [0, 5, 10, 15, 20, 25, 30, 35, 50, 75, 100];%Time
Temp_in = [10, 4, 2, 15, 15, 2, 10, 2, 2, 15, 10]; %Time
dt = 0.1; %stepsize
n = max(t_in)/dt+1; %DataPoints

% parameters
k = 3.76e-3;
gamma = 1.0;
a = 8.86e-2;
Tmin = 8.91;
Ymax = 8.86; %log10(Nmax)
Nmax = 10^Ymax;

f = @(t,N)[-gamma*rate*N(1);gamma*rate*N(1)*m + rate*N(2)*(1-m*(N(1)+N(2))/Nmax)];   %lag cells and dividing cells

% Intial conditions
t(1) = 0;
N(:,1) = [100,0];%intial cell numbers，N1=1000和N2=0
    function umax = rate(i)
        %ti = linspace(0, max(t_in), n);
        Temp = interp1(t_in,Temp_in,n);% T at ti0,1,2..n
        i = 1:n;
           if Temp(i) < Tmin
              umax = k*(Temp(i) -Tmin);
           else
              umax = a^2*(Temp(i) -Tmin)^1.5;
           end
        %end
    end
    
    function m
        i = 1:n;
        if Temp(i) < Tmin
              m = 0;
        else
              m = 1;
        end
    end

%RK4 update loop
for i = 1:n
    t(i+1)=t(i)+dt;
    k1 = f(t(i)     ,N(:,i)        );
    k2 = f(t(i)+dt/2,N(:,i)+k1*dt/2);
    k3 = f(t(i)+dt/2,N(:,i)+k2*dt/2);
    k4 = f(t(i)+dt  ,N(:,i)+k3*dt  );
    N(:,i+1) = N(:,i)+dt/6*(k1+2*k2+2*k3+k4);
end
y = log((N(1,:))+N(2,:))/2.303; %N = N1+N2,transfe to log10

%plot the solution
plot(t,y)
xlablel(Time)
ylablel(Number)
end

[Jan]

Please note, that clear all is completely useless on top of a function, but wastes time only. The brute clearing header "clear all, clc, close all" is cargo-cult-programming (link): It might have look useful in another context, but now it has become a useless piece of magic.

" I know there is something wrong with my code" - then please be so kind and share with the readers, what let you assume this. Providing all useful information is a good strategy, if you want others to help you.

What is the purpose of this part:

    function m
        i = 1:n;
        if Temp(i) < Tmin
              m = 0;
        else
              m = 1;
        end
    end

I cannot guess its intention. Because i is a vector, the command if Temp(i) < Tmin is equivalent to:

    if all(Temp(1:n) < Tmin)

Is this really wanted? It is at least strange to use the name of a nested function as output variable. Calling the function "m" will define teh variable "m", such that the function is shadowed for following calls.

You define "rate" as a function, which expects 1 input argument::

        function umax = rate(i)

But then you call it without inputs:

        ... [-gamma*rate*N(1); ...

You are using anonymous and nested functions. This increases the complexity of the code without any need. I recommend to use standard functions to avoid an unnecessary confusion. As soon as f is defined as standard function, interpolating the temperature will be much easier.

[Shanshui Li]

I am very sorry, I have just begaun to learn matlb,and I tried this code, it doesn't work. Thank you all the same!

[Jan]

You are welcome. This is a forum for solving questions concerning Matlab.

[Gorka Bidegain]

Hi Sanshui,

what do you want to represent by the parameter gamma in equation 1 and 2 (gamma=1)?

Thanks

Gorka

Answer 1

[Jan]

A simpler version without nested and anonymous functions:

function main

% To get the temperature by interpolation:
t_in = [0, 5, 10, 15, 20, 25, 30, 35, 50, 75, 100];  % Time
T_in = [10, 4, 2, 15, 15, 2, 10, 2, 2, 15, 10];      % Temperature

dt   = 0.1; %stepsize
n    = t_in(end) / dt + 1; %DataPoints

t = zeros(1, n+1);  % Pre-allocate the output
N = zeros(2, n+1);

% Intial conditions
t(1)   = 0;
N(:,1) = [100; 0];  %intial cell numbers，N1=1000和N2=0

%RK4 update loop
for k = 1:n
    % Obtain the current temperature:
    % This assumes a fixed temperature during the calculation of a step!!!
    T = interp1(t_in, T_in, k);  % T at ti0,1,2..n
       
    t(k+1)   = t(k) + dt;
    k1       = fcn(t(k),        N(:,k),           T);
    k2       = fcn(t(k) + dt/2, N(:,k) + k1*dt/2, T);
    k3       = fcn(t(k) + dt/2, N(:,k) + k2*dt/2, T);
    k4       = fcn(t(k) + dt,   N(:,k) + k3*dt,   T);
    N(:,k+1) = N(:,k) + dt / 6 * (k1 + 2*k2 + 2*k3 + k4);
end
y = log(sum(N, 1)) / 2.303; %N = N1+N2,transfe to log10

%plot the solution
plot(t,y)
xlablel('Time')     % Not: xlablel(Time)
ylablel('Number')   % Not: ylablel(Number)  Quotes are required
end

function dN = fcn(t, N, T)

% parameters
k     = 3.76e-3;
gamma = 1.0;
a     = 8.86e-2;
Tmin  = 8.91;
Ymax  = 8.86; %log10(Nmax)
Nmax  = 10^Ymax;

if T <= Tmin
   umax = k * (T - Tmin);
   m    = 1;
else
   umax = a^2 * (T - Tmin) ^ 1.5;
   m     = 0;
end

dN = [-gamma * umax * N(1); ...
       gamma * umax * N(1) * m + umax * N(2) * (1 - m * (N(1) + N(2)) / Nmax)];
end

This approach uses a fixed temperature at each step of the solver. You can move the interpolation into the function to be integrated also, to consider the time points exactly, when the temperature is changed immediately:

function main

% To get the temperature by interpolation:
t_in = [0, 5, 10, 15, 20, 25, 30, 35, 50, 75, 100];  % Time
T_in = [10, 4, 2, 15, 15, 2, 10, 2, 2, 15, 10];      % Temperature

dt   = 0.1; %stepsize
n    = t_in(end) / dt + 1; %DataPoints

t_in_k = t_in * dt;

t = zeros(1, n+1);  % Pre-allocate the output
N = zeros(2, n+1);

% Intial conditions
t(1)   = 0;
N(:,1) = [100; 0];  %intial cell numbers，N1=1000和N2=0

%RK4 update loop
for k = 1:n     
    t(k+1)   = t(k) + dt;
    k1       = fcn(t(k),        N(:,k),           t_in_k, T_in);
    k2       = fcn(t(k) + dt/2, N(:,k) + k1*dt/2, t_in_k, T_in);
    k3       = fcn(t(k) + dt/2, N(:,k) + k2*dt/2, t_in_k, T_in);
    k4       = fcn(t(k) + dt,   N(:,k) + k3*dt,   t_in_k, T_in);
    N(:,k+1) = N(:,k) + dt / 6 * (k1 + 2*k2 + 2*k3 + k4);
end
y = log(sum(N, 1)) / 2.303; %N = N1+N2,transfe to log10

%plot the solution
plot(t,y)
xlablel('Time')     % Not: xlablel(Time)
ylablel('Number')   % Not: ylablel(Number)  Quotes are required
end

function dN = fcn(t, N, t_in_k, T_in)

% parameters
k     = 3.76e-3;
gamma = 1.0;
a     = 8.86e-2;
Tmin  = 8.91;
Ymax  = 8.86; %log10(Nmax)
Nmax  = 10^Ymax;

% Obtain the current temperature:
T = interp1(t_in_k, T_in, t);
if T <= Tmin
   umax = k * (T - Tmin);
   m    = 1;
else
   umax = a^2 * (T - Tmin) ^ 1.5;
   m     = 0;
end

dN = [-gamma * umax * N(1); ...
       gamma * umax * N(1) * m + umax * N(2) * (1 - m * (N(1) + N(2)) / Nmax)];
end

[Shanshui Li]

Thank you very much, @Jan!That's the very RK4 method which I was looking for a long time. I am really appreciated your help. I learned much more than the question I asked. By the way, in my question, when T <= Tmin, m equals 0,and when T>=Tmin,m equals 1, but that doesn't matter, I can correct it.

[Gorka Bidegain]

HI, is this model continuous time model solved by rk4? or is it a discrete time model ?

Thanks

[Torsten]

rk4 solves ordinary differential equations in continuous form.