4th order runge kutta for spring mass sytem
5 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Dhrumil Patadia
il 4 Lug 2020
Commentato: Dhrumil Patadia
il 27 Lug 2020
What is wrong with the code: (Also, I am a beginner, so any suggestions for where can i learn matlab for solving odes and pdes without using ode45)
function rk()
Fo=1;m=1;wn=1;w=0.4;z=0.03;
f1=@(t,y1,y2)(y2);
f2=@(t,y1,y2)(Fo/m*sin(w*t)-2*z*wn*y2-wn*y1^2);
h=0.01;
t(1)=0;y1(1)=0;y2(1)=0;
for i=1:10000
t(i+1)=t(i)+h;
k1y1 = h*f1(t(i), y1(i), y2(i));
k1y2 = h*f2(t(i), y1(i), y2(i));
k2y1 = h*f1(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);
k2y2 = h*f2(t(i)+h/2,y1(i)+k1y1/2,y2(i)+k1y2/2);
k3y1 = h*f1(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);
k3y2 = h*f2(t(i)+h/2,y1(i)+k2y1/2,y2(i)+k2y2/2);
k4y1 = h*f1(t(i)+h, y1(i)+k3y1, y2(i)+k3y2);
k4y2 = h*f2(t(i)+h, y1(i)+k3y1, y2(i)+k3y2);
y1(i+1)=y1(i) + (k1y1 + 2*k2y1 + 2*k3y1 + k4y1)/6;
y2(i+1)=y2(i) + (k1y2 + 2*k2y2 + 2*k3y2 + k4y2)/6;
end
plot(t,y1)
5 Commenti
Risposta accettata
Alan Stevens
il 25 Lug 2020
I think your definition of f2 is causingthe problem. Shouldn't it be:
f2=@(t,y1,y2)(Fo/m*sin(w*t)-2*z*wn*y2-wn*y1*abs(y1));
3 Commenti
Alan Stevens
il 25 Lug 2020
When y1 is positive there is no difference; however, when y1 is negative, y1^2 is posiive, but y1*abs(y1) is negative.
Più risposte (0)
Vedere anche
Categorie
Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!