Adapt RK4 ODE Solver from first order to 2nd order
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Nicolas Caro
le 7 Mai 2021
Modifié(e) : Nicolas Caro
le 17 Mai 2021
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 commentaires
James Tursa
le 7 Mai 2021
Modifié(e) : James Tursa
le 7 Mai 2021
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?
Réponse acceptée
James Tursa
le 7 Mai 2021
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 commentaires
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!