Implementation of Runge Kutta Numerical Solution for system of ODE's

조회 수: 8(최근 30일)
Rodrigo Pena 2021년 6월 22일
댓글: Rodrigo Pena 2021년 7월 19일
Hello All,
I am trying to solve a system of ODE's using Runge Kutta method:
I don't know why t, u, x, y and w are not being created ( these are the answers )
Here it is my code:
%Equation 32
% Initial Conditions.
x0 = 0;
t0 = 0;
h = 0.01; % Stepsize.
n = 16; % Number of steps.
function [t,x] = RK21(xdot = u,t0,x0,h,n);
x = [x0];
t = [t0];
for i = 1:n
k1x = f(t(i) , x(i));
k2x = f(t(i) + h, x(i) + k1x*h);
x = [x ; x(i)+ h/2 * (k1x + k2x)];
endfor
endfunction
fi = 45.494 * pi/180; % Longitude.
wz = 7.292*10^-5 * sin(fi); % This is omega z.
g = 9.806; % Gravity [m/s^2]
l = 67; % Lenght [m]
%Equation 34
% Initial Conditions.
u0 = 0;
t0 = 0;
function [t,u] = RK22(udot = 2*wz*w - sqrt(g/l)*x,t0,u0,h,n);
u = [u0];
t = [t0];
for i = 1:n;
k1u = f(t(i) , u(i));
k2u = f(t(i) + h, u(i) + k1u*h);
u = [u ; u(i) + h/2 * (k1u + k2u)];
endfor
endfunction
% Equation 33
% Initial Conditions.
y0 = 0;
t0 = 0;
function [t,y] = RK23 (ydot = w,t0,y0,h,n);
y = [y0];
t = [t0];
for i = 1:n;
k1y = f(t(i) , y(i));
k2y = f(t(i) + h, y(i) + k1y*h);
y = [ y ; y(i) + h/2 * (k1y + k2y)];
endfor
endfunction
% Equation 35
% Initial Conditions.
w0 = 0;
t0 = 0;
function [t,w] = RK24 (wdot = -2*wz*u - sqrt(g/l)*y,t0,w0,h,n);
w = [w0];
t = [t0];
for i = 1:n;
k1w = f(t(i) , w(i));
k2w = f(t(i) + h, w(i) + k1w*h);
w = [ w ; w(i) + h/2 * (k1w + k2w)];
endfor
endfunction
These are the equations that I am trying to solve:
Equations 32, 33, 34 and 35 are first order ODE's from equations 30 and 31. 36, 37, 38 and 39 are the initial conditions that I am using.

댓글을 달려면 로그인하십시오.

채택된 답변

James Tursa 2021년 6월 22일
편집: James Tursa 2021년 6월 22일
The main technical problem (other than your syntax not being MATLAB) is that you are trying to solve your four 1st order equations separately. You can't do this since they depend on each other. You need to solve all four 1st order equations simultaneously. That is, you need to have only one loop and inside that loop you need to use one derivative function that takes a 4-element state vector as input (the x, y, u, w) and outputs a 4-element derivative vector (the dx/dt, dy/dt, du/dt, dw/dt).
You should also pre-allocate your result and then assign elements to this result within your loop. The way you have it currently coded incrementally builds the size of the result within the loop, which causes a lot of unnecessary data copying and can be R E A L L Y S L O W.
댓글 수: 1표시숨기기 없음
Rodrigo Pena 2021년 7월 19일
Thanks !

댓글을 달려면 로그인하십시오.

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by