# Adapt RK4 ODE Solver from first order to 2nd order

조회 수: 1(최근 30일)
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표시숨기기 이전 댓글 수: 2
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 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표시숨기기 이전 댓글 수: 2
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
end

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

### Community Treasure Hunt

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

Start Hunting!