Effacer les filtres
Effacer les filtres

Runge-Kutta 3 variables, 3 equations

20 vues (au cours des 30 derniers jours)
jsparkes951
jsparkes951 le 4 Mai 2015
Commenté : SHIVANI TIWARI le 15 Mai 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 commentaires
Miguel Buitrago
Miguel Buitrago le 13 Juil 2016
Modifié(e) : Miguel Buitrago le 13 Juil 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 le 15 Mai 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

Connectez-vous pour commenter.

Réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by