二阶微分方程龙格库塔法_二阶龙格库塔法公式

二阶微分方程龙格库塔法_二阶龙格库塔法公式一、计算公式对于形如以下的常微分方程:采用四阶龙格库塔法的计算公式:四阶龙格库塔法精度为4,属于单步递推法。单步递推法的基本思想是从(x(i),y(i))点出发,以某一斜率沿直线达到(x(i+1),y(i+1))点。二、实例计算对于下述二阶方程:f(q)为分段函数1、基本思想:令位移q为y(1),q的一阶导dq/dt为y(2),因此可得:dq/dt=y(2)令f(q)=fy令q的二阶导ddq/dt^2=-2*eptheton*y(2)-fy+Fm

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全家桶1年46,售后保障稳定

一、计算公式

对于形如以下的常微分方程:

\frac{dy}{dx}=f(x,y)

y(a)=y0

采用四阶龙格库塔法的计算公式

y_{i+1}=y_i+h/6*(K_1+2K_2+2*K_3+K_4)

K_1=f(x_i,y_i)

K_2=f(x_i+h/2,y_i+h/2*K_1)

K_3=f(x_i+h/2,y_i+h/2*K_2)

K_4=f(x_i+h,y_i+h*K_3)

四阶 龙格库塔法精度为4,属于单步递推法

单步递推法的基本思想是从(x(i),y(i))点出发,以某一斜率沿直线达到(x(i+1),y(i+1))点。

二、实例计算

从上述定义可以看出,龙格库塔实质上是求一阶微分方程,但是如果将一阶导看作变量,则二阶导也不过是这个变量的一阶导而已。

接下来就看实例吧:

对于下述二阶方程:

{\ddot{q}}+2*\varepsilon *{\dot{q}}+q=Fm+Fah*cos(\Omega *t)

1、基本思想:

令位移为    q=u(1)

q的一阶导,即位移的一阶导(速度)为   \dot{q}=u(2)

q的二阶导   \dot{u(2))}=-2* \varepsilon *u(2)-u(1)+Fm+Fah*cos(\Omega *t)

求解位移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);

Jetbrains全家桶1年46,售后保障稳定

求解速度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

运算结果图如下所示

二阶微分方程龙格库塔法_二阶龙格库塔法公式

两中方法计算的结果是一样的。

创作不易,如有帮助,记得点个赞呐。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/234982.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • python可以自动回收垃圾吗_python 数据清洗

    python可以自动回收垃圾吗_python 数据清洗前言现在的高级语言如java,c#等,都采用了垃圾回收机制,而不再像c,c++里,需要用户自己管理内存。自己管理内存及其自由,可以任意申请内存,但这如同一把双刃剑,可能会造成内存泄漏,空指针等bug

  • 几种常见的损失函数「建议收藏」

    几种常见的损失函数「建议收藏」1.损失函数、代价函数与目标函数  损失函数(LossFunction):是定义在单个样本上的,是指一个样本的误差。  代价函数(CostFun

  • GitLab分支介绍

    GitLab分支介绍GitLab简介:GitLab是由GitLabInc.开发,使用MIT许可证的基于网络的Git仓库管理工具,且具有wiki和issue跟踪功能。使用Git作为代码管理工具,并在此基础上搭建起来的web服务。GitLab由乌克兰程序员DmitriyZaporozhets和ValerySizov开发,它由Ruby写成。后来,一些部分用Go语言重写。截止2018年5月,该公司约有290名团队成员,以…

  • java标记符_java标识符的书写规则

    java标记符_java标识符的书写规则标识符就是给java中的类,方法,变量起名字的符号。规则:1.可以由字母,数字,下划线,美元符号组成,但是不能包含@,%,空格等其它特殊字符,不能以数字开头。2.不能是java的关键字,比如static,import,float,int,abstract还有好多这里就不一一例举了,需要的可以自己百度。3.对大小写是敏感的,比如一些包名和类名以及变量名我们虽然大小写不会报错,但是我们需要尽可能的规范…

  • tabnine 激活码【2021免费激活】

    (tabnine 激活码)这是一篇idea技术相关文章,由全栈君为大家提供,主要知识点是关于2021JetBrains全家桶永久激活码的内容https://javaforall.cn/100143.htmlIntelliJ2021最新激活注册码,破解教程可免费永久激活,亲测有效,上面是详细链接哦~DB847YMYYZ-eyJsaWNlb…

  • linux抓包-tcpdump

    linux抓包-tcpdump文章目录1.tcpdump简介2.tcpdump参数3.tcpdump过滤器4.tcpdump常用操作1.tcpdump简介tcpdump是linux平台的抓包工具,可以抓取TCP/IP协议的数据包,网络协议,主机,端口,还提供and,or,not等逻辑语句过滤信息。2.tcpdump参数tcpdump帮助查看tcpdump-h,mantcpdump[root@master~]#tcpdump-htcpdumpversion4.9.2libpcapversion1.5.3

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号