本文介绍如何使用简单的欧拉法求解微分方程，大部分内容出自吴一东老师在他的B站个人空间发布的课程
方法介绍
对于一个一般的微分方程：
{dydt=f(y(t),t)y(0)=y0\begin{cases} \begin{aligned} \frac{\mathrm{d} y}{\mathrm{d} t} &= f(y(t), t)\\ y(0) &= y_0 \\ \end{aligned} \end{cases} ⎩⎨⎧dtdyy(0)=f(y(t),t)=y0
假如我们很难得到他的解析解或者不存在解析解，那我们可以尝试使用欧拉法将微分方程利用泰勒展开构造成一个递推式，通过程序将所有的点求出，进而得到他的数值解。以下说明求法：
区间划分
对于一个微分方程我们只关心函数的一段区域内的数值解。
t∈[0,T]t \in [0, T]t∈[0,T] 将区间划成NNN段，则有Δt=TN\Delta t = \frac{T}{N}Δt=NT
t1=0,t2=Δt,...,tn=(n−1)Δt,...,tN+1=Tt_1 = 0, t_2 = \Delta t, ..., t_{n} = (n-1)\Delta t ,...,t_{N+1} = Tt1=0,t2=Δt,...,tn=(n−1)Δt,...,tN+1=T
泰勒展开
y(t+Δt)=y(t)+dydt(t)Δt+O(Δt2)=y(t)+f(y(t),t)Δt+O(Δt2)\begin{aligned} y(t+\Delta t) &= y(t) + \frac{\mathrm{d} y} {\mathrm{d} t}(t)\Delta t + O(\Delta t^2) \\ &=y(t) + f(y(t), t)\Delta t + O(\Delta t^2) \end{aligned} y(t+Δt)=y(t)+dtdy(t)Δt+O(Δt2)=y(t)+f(y(t),t)Δt+O(Δt2)
得到递推公式
y(tn+1)=y(tn+Δt)=y(tn)+f(y(tn),tn)Δt+O(Δt2)\begin{aligned} y(t_{n+1}) &= y(t_n + \Delta t)\\ &=y(t_n) + f(y(t_n), t_n)\Delta t + O(\Delta t^2) \end{aligned} y(tn+1)=y(tn+Δt)=y(tn)+f(y(tn),tn)Δt+O(Δt2)
近似计算
y(tn)→tnyn+1=yn+f(yn,t)Δty(t_n) \to t_n\\ y_{n + 1} = y_n +f(y_n, t)\Delta t y(tn)→tnyn+1=yn+f(yn,t)Δt
示例
例1：
{f(y,t)=ydydt=yy(0)=1\begin{cases} \begin{aligned} f(y, t) &= y\\ \frac{\mathrm{d} y}{\mathrm{d} t} &= y\\ y(0) &= 1 \end{aligned} \end{cases} ⎩⎨⎧f(y,t)dtdyy(0)=y=y=1
由题目其实我们很容易得知，y=ety = e^ty=et 所以我们可以选择解析解对数值解方法进行验证
按上文求解方法我们求其数值解
yn+1=yn+ynΔty_{n + 1} = y_n + y_n\Delta tyn+1=yn+ynΔt
求解代码如下：
import numpy as npimport matplotlib.pyplot as plt# 设置初始条件T = 5N = 10000dt = T / Nt = np.linspace(0, T, N + 1)# 边界条件y = np.zeros(N + 1)y[0] = 1for i in range(0, N):y[i + 1] = y[i] + y[i] * dtex = np.exp(t)plt.plot(t, y, color='red')plt.plot(t, ex, color='blue')plt.show()
