MATLAB Answers

Solving Lorenz attractor equations using Runge kutta (RK4) method

165 views (last 30 days)
reema shrestha
reema shrestha on 11 Oct 2017
Answered: tyfcwgl on 29 Oct 2017
I am trying to write a code for the simulation of lorenz attractor using rk4 method. Here is the code:
clc;
clear all;
t(1)=0; %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10; %value of constants
rho=28;
beta=(8/3);
h=0.1; %step size
t=0:h:100;
f=@(t,x,y,z) sigma*(y-x); %ode
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
k1=h*f(t(i),x(i),y(i),z(i));
l1=h*g(t(i),x(i),y(i),z(i));
m1=h*p(t(i),x(i),y(i),z(i));
k2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
l2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
m2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
k3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
l3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
m3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
k4=h*f(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
l4=h*g(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
m4=h*p(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
x(i+1)=x(i)+h*(1/6)*(k1+2*k2+2*k3+k4); %final equations
y(i+1)=y(i)+h*(1/6)*(k1+2*k2+2*k3+k4);
z(i+1)=z(i)+h*(1/6)*(m1+2*m2+2*m3+m4);
end
plot3(x,y,z)
But the solutions are not right. I don't know what to do. I know we can do using ode solvers but i wanted to do using rk4 method. I searched for the solutions in different sites but i didn't find many using rk4. While there were some but only algorithm. I tried to compare my solutions with ode45 but doesn't match at all. it's totally different.

  0 Comments

Sign in to comment.

Accepted Answer

Mischa Kim
Mischa Kim on 11 Oct 2017
The only bug that I can see at first glance is here
y(i+1) = y(i) + h*(1/6)*(l1+2*l2+2*l3+l4); % replace ki by li
You also might want to play with (decrease) the step size.

  4 Comments

Show 1 older comment
reema shrestha
reema shrestha on 11 Oct 2017
Thanks a lot. I'd been trying hard and the results were nowhere near to the expected outcome. I've corrected it but still I am not satisfied with the result. It is not same as the result obtained from the ode45.
reema shrestha
reema shrestha on 12 Oct 2017
Wow! Thanks alot. This is beautiful. I wonder why it didn't it worked with mine. Have to learn alot. :)

Sign in to comment.

More Answers (1)

tyfcwgl
tyfcwgl on 29 Oct 2017
I think there are many bugs in your code. After my modifying, it works well. the code and result are below.
clc;
clear all;
t(1)=0; %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10; %value of constants
rho=28;
beta=(8.0/3.0);
h=0.01; %step size
t=0:h:20;
f=@(t,x,y,z) sigma*(y-x); %ode
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
k1=f(t(i),x(i),y(i),z(i));
l1=g(t(i),x(i),y(i),z(i));
m1=p(t(i),x(i),y(i),z(i));
k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
m2=p(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
m3=p(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
k4=f(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
l4=g(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
m4=p(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
x(i+1) = x(i) + h*(k1 +2*k2 +2*k3 +k4)/6; %final equations
y(i+1) = y(i) + h*(l1 +2*l2 +2*l3 +l4)/6;
z(i+1) = z(i) + h*(m1+2*m2 +2*m3 +m4)/6;
end
plot3(x,y,z)

  0 Comments

Sign in to comment.

Sign in to answer this question.

Tags