Adapt RK4 ODE Solver from first order to 2nd order

조회 수: 1 (최근 30일)
Nicolas Caro
Nicolas Caro 2021년 5월 7일
편집: Nicolas Caro 2021년 5월 17일
I have a Matlab function for doing Runge-Kutta4k approximation for first-order ODE's, and I want to adapt it to work for second-order ODE's. This is what I have for first order RK4 ( and it's not working x) ) Would anyone be able to help me find a solution ?
function yt = RK4_2ndOrder(dxdt,y0,h,tfinal)
yt = zeros(2,length(h:tfinal)); %Memory Allocation
yt(:,1) = y0; %Initial condition
i = 2;
for t = h : h : tfinal %RK4 loop
k1 = dxdt(t-h,yt(:,i-1));
k2 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k1*h);
k3 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k1*h);
k4 = dxdt(t, yt(:,i-1) + (k3 * h));
yt(:,i) = yt(:,i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
i = i + 1;
end
end
This is my main
clc;
clear;
h = 0.01; %Time step
y0 = 1; %Initial condition dx1/dt = 1 m.s
tfinal = 20;
tarray = 0:h:tfinal;
ytRK4 = RK4(@dxdt,y0,h,tfinal);
plot(tarray,ytRK4_2ndOrder,'g');
And my function
function dx = dxdt(t,x)
m1 = 12000;
m2 = 10000;
m3 = 8000;
k1 = 3000;
k2 = 2400;
k3 = 1800;
dx = [ x(4) ; x(5) ; x(6) ; -k1/m1*x(1)+k2/m1*(x(2)-x(1)) ; k2/m2*(x(1)-x(2))+k3/m2(x(3)-x(2)) ; k3/m3*(x(2)-x(3)) ];
end
When I try to run my program, it says "Index exceeds the number of array elements" and "Error in dxdt (line 8)
dx = [x(4);x(5);x(6);-k1/m1*x(1)+k2/m1*(x(2)-x(1));k2/m2*(x(1)-x(2))+k3/m2(x(3)-x(2));k3/m3*(x(2)-x(3))];". As a beginner on Matlab, I don't really understand what that means...
Thanks in advance for the help !
  댓글 수: 3
Nicolas Caro
Nicolas Caro 2021년 5월 7일
My bad about RK4 I forgot it :
function yt = RK4(y0,h,tfinal)
yt = y0;
i = 2;
for t = h : h : tfinal
k1 = f(t-h,yt(i-1));
k2 = f(t - (0.5*h), yt(i-1) + 0.5*k1*h);
k3 = f(t - (0.5*h), yt(i-1) + 0.5*k2*h);
k4 = f(t, yt(i-1) + (k3 * h));
yt(i) = yt(i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
i = i + 1;
end
end
I'll call RK4_2ndOrder( ) in the main, like this : (my bad I forgot to copy/paste that line)
ytRK4_2ndOrder = RK4_2ndOrder(@dxdt,y0,h,tfinal);
About yt, I probably messed up because I was trying to understand the logic of other codes to apply it on mine but as you can see, it didn't work ^^.
About dxdt() I'm not sure if I'm doing things correctly, my 3 functions are :
I have to find x1, x2 and x3 and the only initial value I have id dx1/dt (and I don't know how to translate it in my code ^^).
I hope it's a bit clearer now :)
Don't hesitate to tell me if it isn't I'll try my best to explain what I think I'm doing !
James Tursa
James Tursa 2021년 5월 7일
편집: James Tursa 2021년 5월 7일
I can help you set up the code for this, but first you must realize that you HAVE to have the following SIX initial conditions in order to solve this as an initial value problem: x1, x2, x3, x1dot, x2dot, x3dot. You need to recheck your assignment and find all of these initial conditions. Maybe it is just a simple statement that all of the others start at 0?

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

채택된 답변

James Tursa
James Tursa 2021년 5월 7일
Try this for starters
yt = zeros(numel(y0),length(h:tfinal)); %Memory Allocation
And then maybe this is what is intended for initial values:
y0 = [0;0;0;1;0;0]; %Initial condition dx1/dt = 1 m.s
And remember to fix this line for the k2:
k3 = dxdt(t - (0.5*h), yt(:,i-1) + 0.5*k2*h);
  댓글 수: 3
James Tursa
James Tursa 2021년 5월 10일
Please post your current code.
Nicolas Caro
Nicolas Caro 2021년 5월 17일
편집: Nicolas Caro 2021년 5월 17일
function dx = dxdt(t,x)
m1 = 12000;
m2 = 10000;
m3 = 8000;
k1 = 3000;
k2 = 2400;
k3 = 1800;
dx = [x(4);x(5);x(6);-k1/m1*x(1)+k2/m1*(x(2)-x(1));k2/m2*(x(1)-x(2))+k3/m2(x(3)-x(2));k3/m3*(x(2)-x(3))];
end
I didn't see your message earlier sorry, I don't know if i'm doing right to simulate my system with dxdt. For now, what I have is a first order ODE solver using RK4
Main :
clc;
clear;
h = 0.01; %Time step
y0 = 1;
%y0 = [0,0,0,1,0,0]; %Initial condition dx1/dt = 1 m.s
tfinal = 20;
tarray = 0:h:tfinal;
ytRK4 = RK4(y0,h,tfinal);
%ytRK4_2ndOrder = RK4_2ndOrder(y0,h,tfinal);
plot(tarray,ytRK4,'g');
RK4 :
function yt = RK4(y0,h,tfinal)
yt = y0;
i = 2;
for t = h : h : tfinal
k1 = f(t-h , yt(i-1));
k2 = f(t - (0.5*h) , yt(i-1)+0.5*k1*h);
k3 = f(t - (0.5*h) , yt(i-1)+0.5*k2*h);
k4 = f(t , yt(i-1)+(k3 * h));
yt(i) = yt(i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
i = i + 1;
end
end
f
function grad = f(t,y)
grad = (0.1).^t - 3*y;
end

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

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Numeric Solvers에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by