清空筛选条件
清空筛选条件

4th order runge kutta for spring mass sytem

5 次查看（过去 30 天）
Dhrumil Patadia
Dhrumil Patadia 2020-7-4
评论： Dhrumil Patadia 2020-7-27
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 个评论
显示 3更早的评论
Dhrumil Patadia
Dhrumil Patadia 2020-7-25
Why is the number of iteration steps causing a problem. And what if i want solution for a longer time period?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2020-7-25
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 个评论
显示 1更早的评论
Alan Stevens
Alan Stevens 2020-7-25
When y1 is positive there is no difference; however, when y1 is negative, y1^2 is posiive, but y1*abs(y1) is negative.
Dhrumil Patadia
Dhrumil Patadia 2020-7-27
Oh!!
its supposed to be wn^2 * y1, i squared the wrong term. sorry
thanks for your help anyways

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by