Borrar filtros
Borrar filtros

Runge-Kutta 3 variables, 3 equations

20 visualizaciones (últimos 30 días)
jsparkes951
jsparkes951 el 4 de Mayo de 2015
Comentada: SHIVANI TIWARI el 15 de Mayo de 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 comentarios
Miguel Buitrago
Miguel Buitrago el 13 de Jul. de 2016
Editada: Miguel Buitrago el 13 de Jul. de 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 el 15 de Mayo de 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

Iniciar sesión para comentar.

Respuestas (0)

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by