Solving Lorenz attractor equations using Runge kutta (RK4) method

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.

Answer 1

[Mischa Kim]

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.

[Mischa Kim]

reema, I cleaned up the code and vectorized it (speed!).

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.01;   %step size
t=0:h:100;

X      = zeros(3,length(t));
X(:,1) = [x(1); y(1); z(1)];

f = @(t,vec) [sigma*(vec(2) - vec(1));...
              vec(1)*(rho - vec(3)) - vec(2);...
              vec(1)*vec(2) - beta*vec(3)];

for ii = 1:(length(t)-1) %loop
    k1 = f(t(ii)        , X(:,ii)           );
    k2 = f(t(ii) + 0.5*h, X(:,ii) + 0.5*h*k1);
    k3 = f(t(ii) + 0.5*h, X(:,ii) + 0.5*h*k2);
    k4 = f(t(ii) +     h, X(:,ii) +     h*k3);

      X(:,ii+1) = X(:,ii) + h*(1/6)*(k1 + 2*k2 + 2*k3 + k4); %final equations
  end
  plot3(X(1,:),X(2,:),X(3,:))

[reema shrestha]

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.

[Mischa Kim]

See updated answer above.

[reema shrestha]

Wow! Thanks alot. This is beautiful. I wonder why it didn't it worked with mine. Have to learn alot. :)

[Gerardo ramirez]

Can I have your email ? I want to ask you something about attractor.

Answer 2

[tyfcwgl]

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)