引言

考虑存在以下二阶偏微分方程

\[\begin{align}
f_2 \cdot \ddot{X(t)}+f_1 \cdot \dot{X(t)} +f_0 \cdot {X(t)} =F(t)
\end{align}
\]

如何使用四阶龙格-库塔法求解该微分方程?

一阶微分方程的解法

首先回顾下对于一阶微分方程的解法,现在有以下一阶常系数非齐次微分方程

\[\begin{align}
f_1 \cdot \dot{X(t)} +f_0 \cdot {X(t)} =F(t)
\end{align}
\]

方程可以写作

\[\begin{align}

\dot{X(t)} = \frac {F(t)-f_0 \cdot {X(t)}}{f_1 }\\
\frac {X_{n+1}-X_n}{dt} = \frac {F(t)-f_0 \cdot {X(t)}}{f_1 }\\
X_{n+1}=X_n+dt\cdot\frac {F(t)-f_0 \cdot {X_{n}}}{f_1 }=X_n+dt\cdot f(t,X_n)

\end{align}
\]

其中\(dt\)​为时间步长。

四阶龙格-库塔法公式如下:

\[\begin{align}
\begin{cases}

X_{n+1}=X_{n}+\frac{dt}{6}\left (k_1+2k_2+2k_3+k_4 \right) \\
k_1=f\left (t,X_n \right) \\
k_2=f\left (t+\frac {dt}{2},X_n+\frac {dt}{2}\cdot k_1 \right) \\
k_3=f\left (t+\frac {dt}{2},X_n+\frac {dt}{2}\cdot k_2 \right) \\
k_4=f\left (t+{dt},X_n+dt \cdot {k_3} \right)\\

\end{cases}

\end{align}
\]

利用以上公式即可显式推进求解

二阶微分方程的解法

现在考虑二阶常系数非齐次微分方程

\[\begin{align}
f_2 \cdot \ddot{X(t)}+f_1 \cdot \dot{X(t)} +f_0 \cdot {X(t)} =F(t)
\end{align}
\]

方程中出现二次导数项\(\ddot{X(t)}\),难以处理,所以考虑将二次导数项\(\ddot{X(t)}\)转化为一次导数项,至此,我们引入

\[\begin{align}
\begin{cases}

a(t)=X(t) \\
b(t)= \dot a(t)=\dot X(t)

\end{cases}

\end{align}
\]

原方程可以写成一阶偏微分方程组

\[\begin{align}
\begin{cases}

f_2\cdot \dot b(t)+f_1\cdot b(t) +f_0\cdot a(t)= F(t) \\
b(t)=\dot a(t)

\end{cases}

\end{align}
\]

仿照以上一阶微分方程的RK4解法

\[\begin{align}
\begin{cases}
\dot a(t)=b(t)\\
\dot b(t)= \frac {F(t)-(f_1\cdot b(t) +f_0\cdot a(t))}{f_2} \\
\end{cases}
\end{align}
\]
\[\begin{align}
\begin{cases}
\dot a_{n+1}=a_{n}+dt\cdot b(t) = a_{n}+dt\cdot f(t,a_n,b_n) \\

\dot b_{n+1}=b_{n}+dt\cdot \frac {F(t)-(f_1\cdot b(t) +f_0\cdot a(t))}{f_2} = b_{n}+dt\cdot g(t,a_n,b_n) \\
\end{cases}
\end{align}
\]

\[\begin{align}
\begin{cases}
\dot a_{n+1}=a_{n}+dt\cdot f(t,a_n,b_n) \\

\dot b_{n+1} = b_{n}+dt\cdot g(t,a_n,b_n) \\
\end{cases}

\end{align}
\]

化简至此,便可以仿照之前的RK4方法进行求解,具体公式为

\[\begin{bmatrix}
a_{n+1}\\b_{n+1}
\end{bmatrix}=
\begin{bmatrix}
a_n+\frac{dt}{6}\cdot (k_{1a}+2k_{2a}+2k_{3a}+k_{4a})\\
b_n+\frac{dt}{6}\cdot (k_{1b}+2k_{2b}+2k_{3b}+k_{4b})
\end{bmatrix}\\
\]

其中,各个系数求解公式为:

\[\begin{align}
\begin{cases}
\begin{bmatrix}
k_{1a}\\k_{1b}
\end{bmatrix}=
\begin{bmatrix}
f(t,a_n,b_n)\\
g(t,a_n,b_n)
\end{bmatrix}\\

\begin{bmatrix}
k_{2a}\\k_{2b}
\end{bmatrix}=
\begin{bmatrix}
f(t+\frac{dt}{2},a_n+\frac{dt}{2} \cdot k_{1a},b_n+\frac{dt}{2} \cdot k_{1b})\\
g(t+\frac{dt}{2},a_n+\frac{dt}{2} \cdot k_{1a},b_n+\frac{dt}{2} \cdot k_{1b})
\end{bmatrix}\\

\begin{bmatrix}
k_{3a}\\k_{3b}
\end{bmatrix}=
\begin{bmatrix}
f(t+\frac{dt}{2},a_n+\frac{dt}{2} \cdot k_{2a},b_n+\frac{dt}{2} \cdot k_{2b})\\
g(t+\frac{dt}{2},a_n+\frac{dt}{2} \cdot k_{2a},b_n+\frac{dt}{2} \cdot k_{2b})
\end{bmatrix}\\

\begin{bmatrix}
k_{4a}\\k_{4b}
\end{bmatrix}=
\begin{bmatrix}
f(t+dt,a_n+dt \cdot k_{3a},b_n+dt \cdot k_{3b})\\
g(t+dt,a_n+dt \cdot k_{3a},b_n+dt \cdot k_{3b})
\end{bmatrix}\\

\end{cases}

\end{align}
\]

至此已经完成使用四阶龙格-库塔法二阶微分方程求解过程的梳理

实例

求解,以下偏微分方程:

\[\begin{align}
1 \cdot \ddot{X(t)}+0 \cdot \dot{X(t)} +0 \cdot {X(t)} =cos(t)\\
\dot{X(0)}=0\\
{X(0)} =0
\end{align}
\]

理论解:

\[\begin{align}
{X(t)} =-cos(t)+1
\end{align}
\]

代码Gitee仓库链接

#include"RK_ODE.h"
#include<iostream>
using namespace std;
#define M 1
MatrixXd F2(double t) {
MatrixXd res = MatrixXd::Identity(M, M);
return res;
}
MatrixXd F1(double t) {
MatrixXd res = MatrixXd::Zero(M, M);
return res;
}
MatrixXd F0(double t) {
MatrixXd res = MatrixXd::Zero(M, M);
return res;
}
MatrixXd F(double t) {
MatrixXd res = MatrixXd::Ones(M, 1);
res(0, 0) = cos(t);
return res;
}
int main() {
RK_ODE ode(M, 0.1);
ode.X.fill(0.);
ode.dX.fill(0.);
for (int i = 0; i < 1000; i++) {
ode.Solve(F2, F1, F0, F);
ode.WriteMotionAndForce("results1.txt", 0);
}
//ode.SaveSnapshoot("snapshoot.json");
cout << ode._Mode << endl;
}

最新文章

  1. Spark RDD 核心总结
  2. dom4j使用总结
  3. android 使用jdbc1.3.0 操作 sql server
  4. Mysql实现行列转换
  5. excel快捷键设置
  6. Linux进程间通信-匿名管道
  7. 记录C++学习历程
  8. Unity3d之Animation(动画系统)
  9. Qt 学习之路 2(80):定位器
  10. oracle备份表
  11. css内容生成器
  12. Emotional Mastery——英语学习小技巧之一
  13. Eclipse卡顿,内存猛增解决方案
  14. 把VueThink整合到已有ThinkPHP 5.0项目中
  15. C语言链表:删除有序链表中大于mink小于maxk的元素
  16. linux的 .bashrc文件是干什么的?
  17. spring @Validated 注解开发中使用group分组校验
  18. io系列之字符流
  19. 一、HTML基础学习
  20. Vue学习记录(一)

热门文章

  1. LeetCoded第20题题解--有效的括号
  2. ConcurrentModificationException异常原因和解决方法
  3. 事务种类jdbc,Hibernate,JTA事务
  4. Spring系列之集成Druid连接池及监控配置
  5. playwright-python 元素定位、frame处理(一)
  6. 检测一个页面所用的时间的js
  7. tensorflow summary demo with linear-model
  8. AN INTEGER FORMULA FOR FIBONACCI NUMBERS
  9. 揭秘盒马鲜生 Android 短视频秒播优化方案
  10. docker日常使用指南