Solving Lorenz attractor equations using Runge kutta (RK4) method

49 次查看（过去 30 天）
reema shrestha
reema shrestha 2017-10-11
评论： 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.

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

采纳的回答

Mischa Kim
Mischa Kim 2017-10-11
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更早的评论
reema shrestha
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. :)

请先登录，再进行评论。

更多回答（1 个）

tyfcwgl
tyfcwgl 2017-10-29
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)

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by