![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/326789/image.png)
4th order runge kutta for spring mass sytem
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Dhrumil Patadia
am 4 Jul. 2020
Kommentiert: Dhrumil Patadia
am 27 Jul. 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 Kommentare
Akzeptierte Antwort
Alan Stevens
am 25 Jul. 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 Kommentare
Alan Stevens
am 25 Jul. 2020
When y1 is positive there is no difference; however, when y1 is negative, y1^2 is posiive, but y1*abs(y1) is negative.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!