四階龍格-庫塔求二階微分方程

引言

考慮存在以下二階偏微分方程

\[\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;
}