Adapt RK4 ODE Solver from first order to 2nd order
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
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
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 件のコメント
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!