线性常微分方程(Ode)的theta-法(python)

第四十七篇 线性常微分方程的theta-法

线性方程的θ-方法

迄今为止所描述的初值问题数值解的方法，只要它们是一阶导的，并以标准形式排列，适用于线性或非线性方程组。
如果微分方程是线性的，则可能采用的另一种方法为在每一步的开始和结束之间对导数进行线性插值。

一阶方程

考虑一阶线性方程

其中k(x)和l(x)是x的函数。
将x0和x1处的方程微分，距离为h，用简化的符号k(xi) = ki和l(xi) = li，得到

引入了可在此范围内变化的无量纲缩放参数θ

并写出一步法如下

参数θ作为梯度在步骤开始和结束时的权重系数。θ = 0时，为简单欧拉法
最常用的选择是θ = 0.5，它赋予了x0和x1处的梯度相等的权重，等价于积分的梯形法则。在求解时变微分方程组时，使用θ = 0.5有时被称为“Crank-Nicolson”方法。这些θ-方法很受欢迎，因为即使在导数项的额外耦合意味着方程不能容易地简化为方程的“标准形式”也可以使用。
由于方程是线性的，所以很容易得到y1为

这是y1的显式表达以及函数在开始和结束时的h, θ, y0，l0, k0, l1 和 k1 等。对于方程组，分母可以变成一个矩阵，在这种情况下，系统是“隐式的”。

计算实例

h = 0.1使用两步θ-法去计算y(0.2)，θ = 0.5。
由定义的方程得，k(x) =−2 x**2, l(x)=3。
步骤1：
y0 = 0.5, k0 = 0, k1 = −2 × 0.12 = −0.02, l0 = 3, l1 = 3
得到方程

步骤2：
y0 = 0.675, k0 = −0.02, k1 = −2 × 0.22 = −0.08, l0 = 3, l1 = 3
得到方程

这个问题的精确解保留到小数点后三位为y(0.2) = 0.905

二阶方程

利用参数θ进行线性插值也可求解二阶线性方程。
考虑二阶线性方程

将x0和x1处的微分方程，距离h，用简化的符号k(xi) = ki l(xi) = li m(xi) = mi等等，得到

现在得到y1和y’1的下列表达式，利用θ对x0和x1处的导数加权，因此

消去导数项y0 "、y’1和y1 "，重新排列后方程是线性的，因此y1和y’1可成为下面两个方程

从第二个方程的第二项中可以看出，对于θ = 0，该方法将失效。
对于单个二阶方程，四阶龙格-库塔法比本节所述的θ-法要精确得多。然而，对于含有耦合导数方程的大型工程系统，含有θ的线性插值方法由于其简单性，仍被频繁使用

程序如下

求解下面问题