![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/326789/image.png)
4th order runge kutta for spring mass sytem
5 views (last 30 days)
Show older comments
Dhrumil Patadia
on 4 Jul 2020
Commented: Dhrumil Patadia
on 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 Comments
Accepted Answer
Alan Stevens
on 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 Comments
Alan Stevens
on 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.
More Answers (0)
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!