微分方程数值解¶

我们已经知道常系数线性微分方程是比较容易求解的,但是,一旦需要求解非线性的微分方程,我们很多时候都是无能为力,要么是这个微分方程及其难以求解,要么就是这个微分方程干脆没有通解,所以,数值方法就成为我们分析微分方程的一个重要的手段,它不追求求出一个通解,而是采用迭代的手段求出某个值的近似数值解,即我输入一个自变量的值,可以返回给我一个对因变量值的估计值.

欧拉方法¶

考虑一个初值问题的常微分方程:

我们的想法是,采用差分的办法代替微分,不妨设我们要研究的自变量区域为\([a,b]\),把这个区域划分为\(N\)个部分:

步长为:

然后我们用向前差商估计某个点处的导数值:

于是就得到:

我们知道\(y(a)=y_0\),根据这个初值条件迭代可以得到所有的估计值,将\(y(x_n)\)的估计值写做\(y_{n}\),所以就得到近似后的结果:

这样就完成了对微分方程的近似,接下来,一个是可以尝试去求解上述差分方程,另一个则是我们采用迭代的策略得到一系列的离散点,可以直接使用这一系列点,也可以尝试去进行插值和拟合.

如果说采用向后差分,那么近似公式就为:

这个时候我们发现,递推公式不能直接使用了,通常采用所谓试差的方式去计算\(y_{n+1}\),假如已知的是\((x_n,y_n)\),我们先用\(y_n\)估计一下\(y_{n+1}\):

然后把这个估计值继续往递推公式里面去代入,直到结果没有太大的改变

误差估计¶

显然,上述的每一步近似都会带有一定的误差,如何对上面的误差进行一定的估计?

假设点\((x_n,y_n)\)是没有误差的,那么:

定义局部截断误差为:

在\(x=x_n\)处进行泰勒展开:

所以:

所以说,欧拉方法有所谓的二阶截断误差

改进的欧拉方法¶

注意到:

如果使用矩形进行数值积分,就是我们之前使用的欧拉方法,为了减少误差,这次使用梯形进行数值积分,与此同时,估计值的记号和之前一致:

注意到这里的\(y_n\)也不是显式的,所以也要进行迭代试差:

龙格-库塔方法¶

根据拉格朗日中值定理:

不难看出,所谓的近似就是对平均斜率\(y'(\theta h+x_n)\)的估计,普通的欧拉方法使用的是左端点的斜率来估计他,改进了之后的方法就是使用两个端点斜率的平均值来估计他.

所以,龙格-库塔方法的思想就是,在区间内尽可能多取几个点,然后用这几个点的斜率来估计总的平均斜率.

二阶龙格-库塔方法¶

二阶就是取两个点去估计,其中一个点是左端点,另一个就是区间内取的一个点:

其中,可以证明,当这些参数满足下列关系的时候,其阶段误差最小:

4阶龙格-库塔方法¶

4阶龙格库塔方法又称RK方法,比起二阶他多取了两个点,显然精度会更高:

可以调整这些参数使得截断误差尽可能的小,所以这些参数会有限制关系,同理,我们也能够给出一组比较简单的参数满足这个关系,其具体形式如下:

一阶微分方程组:¶

可以通过建立向量关系吧一阶微分方程组的初值问题转化成一阶微分方程的初值问题:

所以一阶微分方程组就转变为:

由于自变量只有x,所以依旧可以采用之前的数值解的方法

高阶微分方程¶

做一个变量代换,把它转变成为一阶微分方程组,令\(y_1=y\):

得到微分方程组:

其初值条件为:

python实现¶

以单摆的运动方程为例,当单摆不满足小角度近似的时候,其方程为:

这是一个非线性的微分方程,我目前只能猜他的解,其解肯定为周期性,但是其周期应该与摆角有关且不初等,所以我们比较难处理,所以采用数值解法,定义角速度:

所以我们就得到微分方程组:

求解器传递给函数的是一个因变量向量和自变量t,函数返回的是一个向量,也就是等式右侧的函数值向量

可以调用库求解:

from scipy.integrate import solve_ivp

总的python代码如下

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

from math import radians
# 常数:
g=9.81
l=1

# 定义微分方程

def F(t,Y):# 注意这里传递参数的时候一定是自变量在前
    theta=Y[0]
    omega=Y[1]

    dth_dt=omega
    dom_dt=-g/l*np.sin(theta)

    return [dth_dt,dom_dt]

y0 = [0, radians(10)]

# 时间区间
t_span = [0, 10]  # 根据实际时间区间设置
t_eval = np.linspace(t_span[0], t_span[1], 1000)

# 求解微分方程
sol = solve_ivp(F, t_span, y0, t_eval=t_eval)

# 绘图
plt.plot(sol.t, sol.y[0], label='theta(t)')
plt.plot(sol.t, sol.y[1], label='omega(t)')
plt.xlabel('Time t')
plt.ylabel('Solution')
plt.legend()
plt.show()

偏微分方程数值解¶

与常微分方程的有限差分法不同的是,偏微分方程存在许多个微分变量,所以不仅仅要在一维数轴上做划分,而是要在二维甚至高维空间上做划分,并且迭代的变量也会增多.

不妨设我们要求的微分方程的解为\(u(x,t)\)

在x轴上做划分,在\(x_j\)处进行泰勒展开:

略去截断误差,并且将函数值改为估计值,就得到向前差分的近似表达式:

以此类推,还有向后差分形式:

以及一阶中心差分:

一维对流方程¶

以一维对流方程为例:

初值条件为:

首先要将空间进行离散化,以便于我们的有限差分的进行:

对x方向和t方向我们采用不同的差分手段,在x方向上采用一阶中心差分,在t方向上采用向前差分:

首先要在区域上进行网格划分:

所以相应的差分形式为:

带入到偏微分方程中得到差分方程:

一维热传导¶

考虑一根一维的铁棒,两端温度维持100℃,在t=0时刻开始向中间进行热传导,试求不同时间下铁棒的温度分布:

一维热传导方程为:

初值条件:

对于二阶导数,我们考虑一阶向后差分为:

所以二阶中心差分为:

对时间维度,采用一阶向前差分:

所以得到递推关系:

整理一下就是:

编写程序来实现这个问题,首先把x和y的取值抽象成为一个个整数,简化问题,例如,t的范围是[0,2.5],x的范围是[0,10],那么就取一个矩阵,大小是25×100.下面我们简化问题,采用20×20的矩阵,并且可以观察不同的导热参数下热传导的变化情况:

import numpy as np
import matplotlib.pyplot as plt


#热传导参数

b=0.5

# 确定步长:

N=20


# 定义一个二维数组,代表函数值

data=np.zeros((N+1,N+1))

for i in range(N+1):
    data[i][N]=100
    data[i][0]=100

for j in range(1,N):
    data[0][j]=0


for n in range(1,N+1):
    for i in range(1,N):
        data[n][i]=data[n-1][i]+b*(data[n-1][i+1]+data[n-1][i-1]-2*data[n-1][i])


# 创建伪彩色图
plt.figure(figsize=(8, 6))
plt.pcolormesh(data, cmap='viridis',shading='auto')
plt.colorbar(label='Value')
plt.title('Pseudocolor Plot')
plt.xlabel('X-axis')
plt.ylabel('Y-axis')
plt.show()

图1.β=0.05	图2.β=0.3	图3.β=0.5

从这几个图可以清晰看到β越大,传热越快,温度梯度的形成更加明显.