![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/326789/image.png)
4th order runge kutta for spring mass sytem
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Dhrumil Patadia
le 4 Juil 2020
Commenté : Dhrumil Patadia
le 27 Juil 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 commentaires
Réponse acceptée
Alan Stevens
le 25 Juil 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 commentaires
Alan Stevens
le 25 Juil 2020
When y1 is positive there is no difference; however, when y1 is negative, y1^2 is posiive, but y1*abs(y1) is negative.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!