【学习总结】动力学方程的龙格库塔积分法（含具体例子与代码）

本文仅用于个人学习总结，如有错误请批评指正。

参考资料

徐超江等，常微分方程基础教程，高等教育出版社，2023年。

1、欧拉法

1.1 前向欧拉

欧拉积分部分不用展开介绍，较为简单。直接拍照课本。

在这里插入图片描述

1.2 梯形法/隐式格式的迭代计算

欧拉法是左侧的数值作为”高度“，所以是一阶的，精度不高。要想高，采用梯形法。

在这里插入图片描述
梯形法是隐式求解，和后面龙格库塔的思路有些相像，所以摆在了这里。

2、龙格库塔积分

先放课本：

在这里插入图片描述

在这里插入图片描述
课本给出了Runge-Kutta的提出思路：避免Taylor展开积分的高阶项造成计算量过大，而是采用这个区间的几个点来对积分进行估计，根据截断误差的阶数，求解对应的系数。可以看出系数是递归计算出来的。

3、高阶常微分方程组的求解

上述例子，都是对于：一元一阶微分方程进行的积分，如果出现了多个变量、或者高阶微分，怎么处理？

处理方法，就是引入中间变量，将高阶微分方程变为多个一阶微分方程，进行求解。参考书籍如下：

在这里插入图片描述注意这里的 $f (x, y, z, w)$ ，就是最高阶数求导对应的方程。仔细观察这组式子，每个等式左侧都是变量/中间变量的一阶导数，右侧都没有导数，因此可以直接积分。

具体的积分方法，就套用上面的Runge-Kutta方法即可，写成下方这种向量形式只是表达简练，但计算时依然每个方程依此计算。

在这里插入图片描述

4、具体例子

下面给出一个具体的动力学方程计算例子：

在这里插入图片描述
可以看出，需要计算的是 $\theta$ 和 $x$ 随时间变化的方程，但原方程中出现了二阶导数。为了敲字方便，导数采用一撇 $^{'}$ 进行标识。

定义变量如下： $y=[\theta', x', \theta, x]$ ，则原方程可以全部换成 $y ()$ 进行表示：
$mL^2 \frac{dy_1}{dt} + mL \frac{dy_2}{dt} + mgL y_3 = 0$ $\frac{dy_1}{dt} + (m+M) \frac{dy_2}{dy} + c y_2 + k y_4 = F(t)$ 以及两个中间变量的定义方程： $\frac{dy_3}{dt}=y_1$ $\frac{dy_4}{dt}=y_2$

这时候，为了凑出来龙格库塔积分时，方程的形式，即左侧是变量的一阶导数，右侧没有微分项。3式和4式是符合的，但1和2并不是。我们发现，如果1式将除了 $\frac{dy_1}{dt}$ 全部外都移动到右侧，右侧存在一个 $\frac{dy_2}{dt}$ ，并不符合基本形式。

因此，我们对1式进行变换，得到： $\frac{dy_1}{dt}= -\frac{1}{L} (\frac{dy_2}{dt} +g y_3)$ ，将这个式子带入2式，得到： $\frac{dy_2}{dt} = \frac{1}{M}(mgy_3-cy_2-ky_4+ \frac{F(t)}{L})$ 同理，利用1式求出 $\frac{dy_2}{dt}=-L\frac{dy_1}{dt}-gy_3$ 后带入2式，得到： $\frac{dy_1}{dt}=\frac{1}{ML}(cy_2+ky_4-(m+M)gy_3- \frac{F(t)}{L})$

此时，修改后的2个式子符合了上述龙格库塔积分的形式。
之后带入书本的式(8.46)，依此计算出第 $n$ 步的 $y_{i,n}$ 即可。

Matlab代码

main函数如下：

clc; clear; close all;
y_init = [0, -1, 0, 0];
tb = 0;
te = 50;
step = 0.1;
[T, R] = rk4(@my_ode, tb, te, y_init, 1000);
figure(1);
plot(T, R(3,:), LineWidth=1)
xlabel("t");
ylabel("\theta");

其中rk4即Runge-Kutta的4阶积分。具体实现如下：

# rk4 的具体计算方法
function [T, R] = rk4(f, a, b, y_init, M)
	h = (b-a)/M;
	Y = zeros(4, M+1);
	T = a:h:b;
	Y(:, 1) = y_init';
	for j = 1:M
		k1 = h*feval(f, T(j), Y(:, j));
		k2 = h*feval(f, T(j) + h/2, Y(:,j)+k1/2);
		k3 = h*feval(f, T(j) + h/2, Y(:,j)+k2/2);
		k4 = h*feval(f, T(j) + h, Y(:, j)+k3);
		Y(:, j+1) = Y(:,j)+(k1 + 2*k2 + 2*k3 + k4)/6;
	end
	R = Y;
end

其中，my_ode 为具体的微分方程：

function eq = my_ode(t, y)
	m=0.4;
	M=1;
	L=0.5;
	g=9.81;
	k=1;
	c=1;
	# da/dt
	dy1_dt = (c*y(2)/M + k*y(4)/M - m*g*y(3)/M - g*y(3) - ft(t)/(M*L)) / L;
	dy2_dt = m*g*y(3)/M - c*y(2)/M - k*y(4)/M + ft(t)/M/L;
	dy3_dt = y(1);
	dy4_dt = y(2);
	eq = [dy1_dt; dy2_dt; dy3_dt; dy4_dt];
end

# ft为外部激励
function y = ft(x)
	y = exp(-x/2);
end