理论力学学习笔记：基于Runge-Kutta法求解非线性微分方程

我们知道，在基于有限元法求解动力学问题时，Newman-β法是最为常用的数值方法。

式中：M为质量矩阵，C为阻尼矩阵，K为刚度矩阵。

Newmark-β法假设在时间间隔[t, t+Δt]内加速度呈线性变化，它的基本假设为：

式中：β和γ为按积分的精度和稳定性要求可以调整的参数。研究表明，当满足

时Newmark-β算法无条件稳定。一般取β=0.5，γ=0.25，此时也称为平均加速度法。

由上述基本假设，可以导出用于迭代的计算式为：

已知t+Δt时刻的振动方程为：

将t+Δt时刻的速度和加速度代入上式，可得：

在实际计算时，一般会先通过第(6)式直接计算得到t+Δt时刻的位移，然后再通过第(4)式中的第一式计算t+Δt时刻的加速度，最后通过第(2)式中的第一式计算t+Δt时刻的速度。

Newmark-β法虽然可以求解非线性振动微分方程，但微分方程组必须满足式(1)所示的形式，对于一些不满足式(1)所示形式的非线性微分方程，则不容易通过Newmark-β法进行求解。此时，可采用Runge-Kutta法进行求解。

Runge-Kutta法是一类显式的单步迭代方法，在实际应用中，一般会采用具有较高精度的四阶Runge-Kutta法，考虑如下所示的一阶常微分方程的初值问题：

其四阶Runge-Kutta法的迭代格式为：

上述方法只能求解一阶常微分方程，对于本文讨论的二阶常微分方程，或者更高阶的常微分方程，需要引入新变量将高阶常微分方程的初值问题转换为一阶常微分方程组的初值问题。考虑如下所示的二阶常微分方程的初值问题：

引入新变量v，并令：

则原二阶常微分方程变为：

上式为一阶常微分方程组，在形式上与单个一阶常微分方程并无不同，只需将该问题记为向量形式：

则原方程组可以写为：

将四阶Runge-Kutta法中的函数换为向量，可得：

如果记

则针对二阶常微分方程的迭代格式的分量形式为：

以上即为针对二阶常微分方程的四阶Runge-Kutta法迭代格式。从上式可以看出，Runge-Kutta法为单步算法，即当前时刻的预测值只需基于上一时刻的预测值进行计算。利用时间t_n上的值u_n和v_n，顺序计算L₁，M₁，L₂，M₂，L₃，M₃，L₄和M₄，即可计算下一时刻t_n+1上的u_n+1和v_n+1。

对于更高阶的常微分方程，其处理方法与上述步骤完全相同。

下面以一道动力学题目为例介绍如何通过Runge-Kutta法求解系统的响应。

均质半圆柱的质量为m，半径为r，可在水平面上作纯滚动。开始时半圆柱的直径AB在铅锤位置，然后半圆柱无初速地向右侧滚动。试建立半圆柱做纯滚动的运动微分方程。

由于半圆柱做纯滚动，因此该系统为单自由度系统。首先需要计算均质半圆柱的质心坐标和转动惯量。

假设半圆柱处于水平位置时的质心坐标为(xC,yC)，如下图所示。根据对称性可得半圆柱在x方向的质心坐标为：

设均质圆柱的面密度为ρ，在该均质圆柱上取一厚度为dy的微段，则该微段的质量为：

已知该微段在y方向的质心坐标为y，因此均质半圆柱在y方向上的质心坐标可表示为：

由于半圆柱绕质心C转动的转动惯量较难计算，因此不妨先计算均质半圆柱绕圆心O点转动的转动惯量。

已知转动惯量的定义为：

因此，以半圆柱的圆心做原点建立极坐标系，并在均质半圆柱上选取一圆弧微段，如下图所示。

已知该圆弧微段的质量可表示为：

因此该圆弧微段对O点的转动惯量为：

因此半圆柱对O点的转动惯量为：

已知半圆柱的面密度为：

代入可得：

当然，半圆柱对O点的转动惯量也可以看作是完整圆柱对O点的转动惯量的一半，注意完整圆柱的质量应为2m，则半圆柱对O点的转动惯量为：

根据平行移轴定理，可得半圆柱对O点和对质心C点的转动惯量之间满足如下关系：

因此可得半圆柱对质心C点的转动惯量为：

（1）动量矩定理

此题采用对接触点D的动量矩定理进行求解最为简单。

考虑均质半圆轮运动如下图所示的一般位置，此时半圆轮的转角为θ。

由于地面对均质半圆轮的约束作用，因此半圆轮的圆心O将始终在水平线上移动，假定D点为车轮滚动转角θ之后与地面的接触点，不难发现OʹD始终垂直于地面。

根据对接触点D的动量矩定理可得：

因此需要计算半圆柱绕D点的转动惯量。

由上图可得∠BʹOʹD=θ，因此有：

因此根据余弦定理可得：

并且已知：

代入可得：

根据平行移轴定理可得半圆轮对当前速度瞬心D的转动惯量为：

代入可得系统的运动微分方程为：

（2）Lagrange方程

此题也可采用分析力学中的Lagrange方程进行求解。

已知半圆轮在当前瞬时的动能可表示为：

取半圆轮直径与地面垂直为零势能点，过Cʹ向OʹD做垂线并交于E点，则OʹE的长度可表示为：

因此半圆轮的势能可表示为：

对于机械能守恒的保守系统，其Lagrange方程可表示为：

式中：q_j为第j个广义坐标，L为Lagrange函数，可定义为：

因此，该系统的Lagrange函数可写为：

因此有：

因此，根据Lagrange方程可得：

为采用Runge-Kutta法求解半圆柱做纯滚动时的运动微分方程，将其整理为如下形式：

令ω=dθ/dt，则上述微分方程可转换为：

基于上式即可根据第2节中给出的Runge-Kutta法的迭代格式进行求解，求解程序可参见附录。

取半圆轮质量m=50 kg，半径r=1 m，重力加速度g=9.81 m/s²，并取初始时刻圆轮的转角为θ₀=0，初始角速度为ω₀= 0，计算时间t=10 s，时间步长h=0.01 s，基于四阶Runge-Kutta法可以计算得到半圆轮做纯滚动时转角和角速度随时间的变化如下图所示。

从图中可以看到，由四阶Runge-Kutta法给出的计算结果与通过MATLAB中的ODE45函数（可变步长的四阶Runge-Kutta法）得到的计算结果完全吻合，但当计算时间较长时，由MATLAB中的ODE45函数得到的计算结果似乎存在较大的误差。例如，下图给出了时间步长不变，并且取计算时间t=100 s的相图（转角与角速度构成的曲线）。

由于半圆柱做周期运动，因此其相轨迹曲线应该构成一条封闭曲线，这与图中采用四阶Runge-Kutta法给出的计算结果一致，而采用MATLAB中的ODE45算法计算得到的相轨迹曲线表现为随时间变化而逐渐膨胀，表明数值误差在逐渐累积。

附录：四阶Runge-Kutta法求解程序

clear,clc
%程序基于四阶Runge-Kutta法求解二阶非线性微分方程
%****************************模型参数设置**********************************
t0=0;                            %计算初始时间
tf=100;                           %计算终止时间

h=0.01;                             %时间增量

m=50;                            %半圆轮质量
r=1;                             %半圆轮半径
Acc_g=9.81;                      %重力加速度

Theta0=0;                        %初始角度
Omega0=0;                        %初始角速度
%**************************************************************************

%*****************************初始化模型***********************************
%计算增量步数
Num_inc=(tf-t0)/h+1;
%初始化时间、转角、角速度
t=zeros(Num_inc,1);
Theta=zeros(Num_inc,1);
Omega=zeros(Num_inc,1);
Theta(1)=Theta0;
Omega(1)=Omega0;
%**************************************************************************

%*******************基于四阶Runge-Kutta法求解微分方程**********************
for i=1:Num_inc-1
  %累加时间增量
  t(i+1)=t(i)+h;
  %计算系数
  L1=f(t(i),Omega(i),Theta(i),Acc_g,r);
  M1=g(t(i),Omega(i),Theta(i));
  L2=f(t(i)+h/2,Omega(i)+h/2*L1,Theta(i)+h/2*M1,Acc_g,r);
  M2=g(t(i)+h/2,Omega(i)+h/2*L1,Theta(i)+h/2*M1);
  L3=f(t(i)+h/2,Omega(i)+h/2*L2,Theta(i)+h/2*M2,Acc_g,r);
  M3=g(t(i)+h/2,Omega(i)+h/2*L2,Theta(i)+h/2*M2);
  L4=f(t(i)+h,Omega(i)+h*L3,Theta(i)+h*M3,Acc_g,r);
  M4=g(t(i)+h,Omega(i)+h*L3,Theta(i)+h*M3);
  %计算下一时刻的转角和角速度
  Omega(i+1)=Omega(i)+h/6*(L1+2*L2+2*L3+L4);
  Theta(i+1)=Theta(i)+h/6*(M1+2*M2+2*M3+M4);
end
%**************************************************************************

%************************绘制转角与角速度随时间变化************************
%绘制转角随时间变化规律
plot(t,Theta);hold on;scatter(t,Theta);grid on;
title("半圆轮转角随时间变化规律");xlabel("时间/s");ylabel("转角/rad");
%绘制角速度随时间变化规律
figure;plot(t,Omega);hold on;scatter(t,Omega);grid on;
title("半圆轮角速度随时间变化规律");xlabel("时间/s");ylabel("角速度/rad/s");
%绘制半圆轮滚动的相图(角度与角速度之间的关系)
figure;plot(Theta,Omega);axis equal;
title("半圆轮做纯滚动的相图");xlabel("转角");ylabel("角速度");
%**************************************************************************

function f=f(t,Omega,Theta,Acc_g,r)
f=8*Acc_g*cos(Theta)/(r*(9*pi-16*sin(Theta)));
end

function g=g(t,Omega,Theta)
g=Omega;
end