runge kutta algae flower
1 view (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Torsten
on 16 Nov 2023
Edited: Torsten
on 16 Nov 2023
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.
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!