显示欧拉法求解 y' = y 问题

简介

最近 无聊看了一下 数值解法

已知 \(\frac{dy}{dx} = y\) 我们知道其有一个解析解为 \(y = e^x\)
同时我们知道 其初值 \(y(0) = 1\), x 的范围为 \(0<= x <= 1\)
我们需要求解其数值解，也就是\(y\) 值在 \(0<= x <= 1\) 的分布

我们猜想显示欧拉是由 泰勒展开式进行构建的

\[f(x) = f(x_0) +f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+···+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x) \]

也就是可以简单等同于求解

\[y(x) = y(0) + \frac{dy}{dx(x == 0)}(x - 0) + 余项（即误差） \]

但是关于 \(\frac{dy}{dx(x == 0)}\) 比较难以求解，采用 \(\frac{\mathbf{y}\left( k+1 \right)-\mathbf{y}\left( k \right)}{h}\) 来近似求解
其中 h 即步长步长越小带来的精度越来，同时带来的计算量越大
这里比较特殊由于题目中已知了\(\frac{dy}{dx} = y\)所以我们可以使用

\[y_{n+1} = y_{n} + h(y') \]

\[y_{n+1} = y_{n} + h y_{n} \]

那么我们对于其数值解就可以比较简单求解出来了

code

import math
from matplotlib import pyplot as plt

t_0 = 0
y_0 = 1
tau = 0.01
i = 1
solve = []
Euler = []
t = []
while i < 1000:
    if i == 1:
        y_n = y_0
        t_n = t_0
    Euler.append(y_n)
    solve.append(math.exp(t_n)) # 精确解
    t.append(t_n)
    func = y_n
    y_n = y_n + tau * func
    t_n = t_n + tau
    i += 1

plt.plot(t, Euler, c='green', label=' Euler method')
plt.plot(t, solve, c='red', label=' accuracy')
plt.fill_between(t, solve, Euler, facecolor='blue', alpha=0.2)
plt.title('Euler method', fontsize=19)
plt.xlabel('t', fontsize=19)
plt.ylabel('y', fontsize=19)
plt.legend()
plt.show()

img

其中紫色部分是误差

参考链接

https://zhuanlan.zhihu.com/p/131303043
https://zhuanlan.zhihu.com/p/422074759