Solving Lorenz attractor equations using Runge kutta (RK4) method

Question

reema shrestha 2017-10-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360751-solving-lorenz-attractor-equations-using-runge-kutta-rk4-method

评论： Gerardo ramirez 2020-5-18

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mischa Kim 2017-10-11

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360751-solving-lorenz-attractor-equations-using-runge-kutta-rk4-method#answer_285234

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

reema shrestha 2017-10-12

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

Gerardo ramirez 2020-5-18

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

请先登录，再进行评论。

Answer 2

tyfcwgl 2017-10-29

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360751-solving-lorenz-attractor-equations-using-runge-kutta-rk4-method#answer_288399

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)