MENU

Numerical Analysis: 微分方程的欧拉法

• February 15, 2020 • Read: 4706 • Knowledge,Python,Numerical Analysis

一阶常微分方程的欧拉法

一阶ODE是最简单的一种ODE,微分方程本身是一个关于函数的方程,如果说这个ODE是显式给出的,那么这实际上给出的就是这个函数的导数的表达式,从数值计算的角度来看,数值方法关心的是整个区域内各个采样点(离散点)处的函数值,而此时显式是函数的导数,那么如果说给这个式子赋值进行计算的话,得到的将是每个采样点处函数的斜率值,这个场也称为斜率场。从斜率场来理解ODE以及Euler法是一个非常直观的。

假设ODE表达式为$y'=ty + t^3$,这个式子是显式表达式,其是一个ODE,同时它也是一个函数导数的表达式,由于方程右侧中不仅含有函数值$y$还包含有自变量$t$,所以这个ODE是非自相关的,同理,如果导数表达式中不含有自变量,就称为自相关的。斜率场就是对这个显式进行可视化的结果,它反映的是各个点处斜率的大小和方向。假设$t$的区间是$[0,1]$,进行斜率场的可视化,这里对于$y$的范围为了方便取$[0,1]$,即结果是反映在$[0,1]$区间上的斜率场分布。

import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0, 1, 11)  # 步长0.1
y = np.linspace(0, 1, 11)
tt, yy = np.meshgrid(t, y)
dy = tt * yy + tt**3
E_t, E_y = np.gradient(dy)
E_norm = np.hypot(E_t, E_y)
fig, ax = plt.subplots()
ax.quiver(tt, yy, E_t/E_norm, E_y/E_norm, color='tab:red')
ax.set_title($''$)
fig.savefig('slope.png', dpi=600)
plt.show()

slopefield

斜率场是认识、研究ODE的一个强有力的工具,也是理解Euler法的基础。试想如果ODE是显式表达的,那么就可以确定出斜率场,它意味着场中各个点处的斜率是已知的,那么对于某一点来说,它的函数值如果已知,然后利用它的斜率,可以计算出函数值变化的增量,进而可以得到下一个点,即这一思路在斜率场中表现为从一点处沿着斜率方向移动到下一个,反复的进行这一过程,即可得到问题的解答,当然这一过程完成,是建立在给定初始值的情况下,也即ODE问题实际上就是一个IVP问题。Euler法就是利用斜率与步长的乘积来确定$y$的增量,利用迭代关系完成整个过程的计算。

$$ w_{0} = y_0\\ w_{i + 1} = w_i + h \times df $$

其中$h$是自变量的步长,即$\mathrm{d}t$,因为ODE是显式表达的,所以$df$是ODE本身。

之前已经说明了Euler法计算的思路,以及斜率场在ODE的IVP问题中的角色。以前述给出的例子进行计算来验证。若设定初始值$y(0) = 0.6$,自变量区间仍是$[0,1]$,进行Euler法的数值实现。

h = 0.1
tmin = 0
tmax = 1
t = np.arange(tmin, tmax+h, h)
w = np.zeros(t.shape[0])  # 初始化一个结果向量
w[0] = 0.6
for i in np.arange(w.shape[0] - 1):
      df = t[i] * w[i] + t[i]**3  # 计算当前点的斜率值
      w[i + 1] = w[i] + h * df  # 计算下一个点的函数值
print(w)

为了体现出Euler法求解IVP问题的核心思想,对于上述方程给定另外一个初始值作为对比、参照,$y(0)=0.4$。结果如下。

slopefield

从斜率场的角度来看,当ODE是一个显式表达时,意味着场中各个离散点处的斜率时已知的,利用函数微分的概念获得函数值的增量,由此得到下一个点处的函数值。而要想完成迭代计算,那么必须提供初始的函数值$y_0$,不同$y_0$意味着从不同的点出发开始移动。

Last Modified: February 16, 2020
Archives Tip
QR Code for this page
Tipping QR Code