Runge-Kutta 3 variables, 3 equations

83 views (last 30 days)
jsparkes951
jsparkes951 on 4 May 2015
Commented: SHIVANI TIWARI on 15 May 2019
I am trying to use the 4th order Runge Kutta method to solve the Lorenz equations over a perios 0<=t<=250 seconds. I am able to solve when there are two equations involved but I do not know what do to for the third one. I have :
x(1)=1;
y(1)=2;
z(1)=0;
h=0.1;
sigma=10;
rho=28;
beta=(8/3);
f=sigma*(y-x);
g=rho*x-y-x*z;
w=x*y-beta*z;
for i=1:2500
k1=h*f(x(i),y(i),z(i));
l1=h*g(x(i),y(i),z(i));
m1=h*w(x(i),y(i),z(i));
I do not know what to do when calculating k2 etc. I have:
k2=f(x(i)+h/2, y(i)+k1/2, z(i)+l1/2);
etc. but is there any change with m(i)? I only know how to solve this for two equations. Am I on the right track? I'm new to MATLAB and Runge-Kutta so any help would be greatly appreciated.
  2 Comments
Miguel Buitrago
Miguel Buitrago on 13 Jul 2016
Edited: Miguel Buitrago on 13 Jul 2016
I used the 4th order Runge Kutta method to solve the Lorenz equations:
k1 = h*f(u1,u2);
m1 = h*g(u1,u2,u3);
n1 = h*e(u1,u2,u3);
k2 = h*f(u1+k1/2, u2+m1/2);
m2 = h*g(u1+k1/2, u2+m1/2, u3+n1/2);
n2 = h*e(u1+k1/2, u2+m1/2, u3+n1/2);
k3 = h*f(u1+k2/2, u2+m2/2);
m3 = h*g(u1+k2/2, u2+m2/2, u3+n2/2);
n3 = h*e(u1+k2/2, u2+m2/2, u3+n2/2);
k4 = h*f(u1+k3, u2+m3);
m4 = h*g(u1+k3, u2+m3, u3+n3);
n4 = h*e(u1+k3, u2+m3, u3+n3);
u1 = u1 + (k1+(k2*2)+(k3*2)+k4)/6;
u2 = u2 + (m1+(m2*2)+(m3*2)+m4)/6;
u3 = u3 + (n1+(n2*2)+(n3*2)+n4)/6;
SHIVANI TIWARI
SHIVANI TIWARI on 15 May 2019
Runge-kutta third order method:
%rk3:runge kutta of thirdorder
clc;
clear all;
close all;
% y' = y-x ode condition
f = @(x,y) y-x;
fex = @(x) exp(x)+x+1; % exact solution
a=0;
b= 3.2;
n =16;
h=(b-a)/n;
y(1) =2; %initial value
i = 0;
for x= a:h:b
i = i+1;
K1 = f(x,y(i)); %initializing solution
K2 = f(x+h*0.5,y(i)+h*K1*0.5);
K3 = f(x+h, y(i)-h*K1 +2*K2*h);
y(i+1) =y(i)+h*(1/6)*(K1 +4*K2+K3);
g(i) = fex(x);
xx(i) = x;
Error(i) = abs(g(i) - y(i)); %error obtain
end
%plot result
plot(xx,y(1:n+1),'k',xx,g,'y')
legend('RK3','Exact solution')
xlabel('x')
ylabel('y')
title('RK3 vs exact solution')
can you plz check my code.actually i dont get exact error. so plz check my prgram

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!