runge kutta algae flower

Hello,

I need help with a home asignment I have for a matlab course. The aim is to solve P and Z for the two functions dP/dt and dZ/dt. I have managed to write the code in ah_alg.m, the function in lorenz and the rk4 in rk4 files. However, I get the wrong answer and as such, my error function is wrong.

The asignment is to a) solve for P and Z. b) Plot the error function for P and Z by calculating |P(n)h/2 - P(n)h| where P(n)h is P with regular steplength, P(n)h/2 is steplength halved.

dP/dt = rP (1 - P2/α2 + P2) - RmZ

dZ/dt = γRmZ P/(K + P) - μZ

This is a breif summary of the asignment, the whole thing can otherwise be found in proj1mat_eng.pdf. As mentioned, my code runs, but produces the wrong result. Can anyone look through my code files and see whats wrong or is it better to start over?

close all 
clc

h = 1;                  % Step length
h2 = h/2;
tspan = 0:h:400;       %Start value, time expressed as days

% Initial conditions
p0 = 20;
z0 = 5;
x0 = [p0,z0]'; % Write everything in vector form
r = 0.3;
K = 108;
R = 0.7;


%Parameters
alpha = 5.7;
mu = 0.024;
gamma = 0.05;

% Solving Runge knutta 4

X(:,1) = x0;
xin = x0;
xin2 = x0;
Err = [0; 0]; % Matrix for saving error for p and x

time = tspan(1);

for i=1:tspan(end)/h
    time = time+h;
    pout = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h,time,xin);
    pout2 = rk4(@(t,x)lorenz(t,x,K,R,r,alpha,mu,gamma),h2,time,xin2);
    X = [X pout];
    Err(1,i) = abs(pout2(1)-pout(1));  %Error for p
    Err(2,i) = abs(pout2(2)-pout(1)); %Error for z
    xin = pout;
    xin2 = pout2;
end


Err = [0 Err(1,:);
       0 Err(2,:)]; % Adding start error, since both functions start at p0,z0 = 0
plot(log(tspan(1,:)),log(Err(1,:))) % Erf for p
hold on
plot(log(tspan(1,:)),log(Err(2,:)))

function dp = lorenz(t,x,K,R,r,alpha,mu,gamma)
  dp = [
  r*x(1)*(1-x(1)/K) - R*x(2)*(x(1)^2/(alpha^2 + x(1)^2));
  gamma*R*x(2)*(x(1)/(alpha^2 + x(1)^2)) - mu*x(2);
  ];
end

function pout = rk4(f,h,tk,xk)
%Rungeknutta 4 method
%   This function calculates Rk4
% Startvalue for k1, t = t0

  k1 = f(tk,xk);  
  k2 = f(tk+h/2,xk+h*k1/2);
  k3 = f(tk+h/2,xk+h*k2/2);
  k4 = f(tk+h,xk+h*k3);
  pout = tk + (h/6)*(k1+k2+k3+k4);
end

Answer 1

[Torsten]

Don't you have to call rk4 twice for stepsize h/2 to compare the results ?

And shouldn't the Runge-Kutta step be

pout = xk + (h/6)*(k1+2*k2+2*k3+k4);

instead of

pout = tk + (h/6)*(k1+k2+k3+k4);

?

And shouldn't

Err(2,i) = abs(pout2(2)-pout(2)); %Error for z

instead of

Err(2,i) = abs(pout2(2)-pout(1)); %Error for z

And shouldn't

time = time+h;

appear at the end of the loop instead of the beginning ?

After these changes, you get reasonable results.

[Alexander]

Thanks for the reply,

I was able to solve for p and z now