MATLAB: How to solve ODEs with time dependent parameters by RK4 method

numerical solutionrk4time dependent odes

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
ODEs.png
Temp-time profile.png

Best Answer

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
Related Question