大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。
Jetbrains全家桶1年46，售后保障稳定
一、计算公式
对于形如以下的常微分方程：
采用四阶龙格库塔法的计算公式：
四阶 龙格库塔法精度为4，属于单步递推法。
单步递推法的基本思想是从（x(i),y(i)）点出发，以某一斜率沿直线达到（x(i+1),y(i+1)）点。
二、实例计算
从上述定义可以看出，龙格库塔实质上是求一阶微分方程，但是如果将一阶导看作变量，则二阶导也不过是这个变量的一阶导而已。
接下来就看实例吧：
对于下述二阶方程：
1、基本思想：
令位移为
q的一阶导，即位移的一阶导（速度）为
q的二阶导
求解位移u(1)的龙格库塔计算方法如下：
KK1=u(2);
KK2=u(2)+h/2*KK1;
KK3=u(2)+h/2*KK2;
KK4=u(2)+h*KK3;
u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);求解速度u(2)的龙格库塔计算方法如下：
K1=-2*eptheton*u(2)-u(1)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K2=-2*eptheton*(u(2)+h/2*K1)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K3=-2*eptheton*(u(2)+h/2*K2)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
K4=-2*eptheton*(u(2)+h*K3)-(u(1)+h)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);
u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);2、编程实现
参数设置：
h=0.001; % 步长
t0=0:h:400; % 总时长
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
u(1)=0;u(2)=0; % 初值总的程序实现
h=0.001;
t0=0:h:400;
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
u(1)=0;u(2)=0;
for i=1:length(t0) % 进行多次迭代
tao=t0(i);
u=RK(u,tao,h,ep,w,Fm,Fah);
Result(i,:)=u; % 将每次迭代的位移和速度保存
end
figure(1)
subplot(2,1,1)
plot(t0,Result(:,1)) % 绘制位移图
xlabel('Time')
ylabel('displacement')
subplot(2,1,2)
plot(t0,Result(:,2)) % 绘制速度图
xlabel('Time')
ylabel('velocity')
function u=RK(u,tao,h,ep,w,Fm,Fah)
KK1=u(2);
KK2=u(2)+h/2*KK1;
KK3=u(2)+h/2*KK2;
KK4=u(2)+h*KK3;
u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);
K1=-2*ep*u(2)-u(1)+Fm+Fah*cos(w*tao);
K2=-2*ep*(u(2)+h/2*K1)-u(1)-h/2+Fm+Fah*cos(w*tao);
K3=-2*ep*(u(2)+h/2*K2)-u(1)-h/2+Fm+Fah*cos(w*tao);
K4=-2*ep*(u(2)+h*K3)-u(1)-h+Fm+Fah*cos(w*tao);
u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);
end结果图如下所示
值得特别注意的是
u=RK(u,tao,h,ep,w,Fm,Fah);输入的u与输出的u一定要符号一致，从而确保下一次迭代的初始值是上一次的值。确保迭代能一直进行下去。
三、辅助验证
接下来用MATLAB自带的ode45函数来进行验证。
之前已经写过ode45函数的用法，在此不再进行介绍。
源程序如下：
t0=0:0.001:400;
w=5;
ep=0.02;
Fm=0.1;
Fah=0.05;
y0=[0 0];
[t,u]=ode45(@(t,u) odefun(t,u,w,ep,Fm,Fah),t0,y0);
figure(1)
subplot(2,1,1)
plot(t,u(:,1))
xlabel('Time')
ylabel('displacement')
subplot(2,1,2)
plot(t,u(:,2))
xlabel('Time')
ylabel('velocity')
function du=odefun(t,u,w,ep,Fm,Fah)
du=[u(2);
-2*ep*u(2)-u(1)+Fm+Fah*cos(w*t)];
end运算结果图如下所示
两中方法计算的结果是一样的。
创作不易，如有帮助，记得点个赞呐。
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/234982.html原文链接：https://javaforall.net
