经典的Euler法提供了认识、求解ODE基本方法和思路，其他的求解方法均是建立在该方法基础上的一个延伸，认识、理解Euler法是非常重要的。在经典Euler法中，函数值的微分$\mathrm{d}y$是依据全微分理论构建得到的，因为ODE自身提供了求解区域中某一点处的斜率值，只要规定好步长，那么对应的函数值增量就得以确定，进而可以得到下一个点处的函数值，反复进行这一过程即可完成求解。而搜寻下一个点的关键就是斜率值，Euler中斜率值是取的当前的点，即$t_i,y_i$处的斜率值，来表示$y_i与y_{i+1}$之间的斜率，如果可以缩小步长，那么精度会上升，另一个角度来看，如果利用该区段中点处的斜率值，可能效果要更好，更具有代表性，基于这一个改进思路就得到改进的Euler法，也称为显示的梯形方法。

Euler法修正的思路

Euler法改进的缘由

Euler法提供的思路是基于斜率场来完成的，对于某一点来说，$df = f(t_i,y_i)'$表示斜率值，那么依据全微分理论，自然可以得到对应区间上的函数值$y$的微分$\mathrm{d}y$值，它反映了$t_i, y_i$到$t_{i+1}, y_{i+1}$的函数值的变化大小，那么由此得到$y_{i+1}$的值。一般来说，$y$值表示ODE的真实值，而数值计算的值是一个估计值（近似值）,为了区别两者，设$w$表示函数值的估计值。

$$ w_{i + 1} = w_{i} + h \times df(t_{i}, w_{i}) $$

前述已经说明了，改进的一种方式就是对斜率值进行修正，经典法中是利用区段左端点的斜率值来作为衡量的标准，如果区段长度过大，精度下降的程度较大。

斜率值的确定

取区间两个端点处斜率的平均值作为该区段的斜率值，并用于计算。

$$ \begin{equation} slope = \frac{1}{2}(df(t_i, w_i) + df(t_{i + 1}, w_{i+1})\\ slope = \frac{1}{2}(df(t_i, w_i) + df(t_{i} + h, w_i + h \cdot df(t_i,w_i))) \end{equation} $$

改进的Euler法，依然延续了经典法的计算思路，由于经典Euler法可以计算得到下一个点处的斜率值，然后取钱一个点和后一个点处的斜率值的平均值作为斜率的修正值，再利用Euler法进行修正后的计算，这就是改进的Euler法，因为该公式与梯形面积的计算公式相似，因此也叫作显式梯形法。

$$ \begin{equation} w_{0} = y_0\\ w_{i+1} = w_{i} + \frac{h}{2} \cdot (df(t_i, w_i) + df(t_i + h, w_i + hdf(t_i, w_i))) \end{equation} $$

通过上述的说明，可以看出改进Euler法本质还是经典的Euler法，只不过对于一次计算，它是在Euler法基础上进行了第二次迭代计算用于修正。

计算验证与对比

依然以之前经典Euler法中的例子进行计算，这样方便对比两者的结果差异。

$$ \begin{equation} y' = t y + t ^ 3 \quad (0\leq t \geq 1)\\ y_0 = 1 \end{equation} $$

对于$t$这里共划分10段，即步长$h=0.1$进行ODE的计算。改进Euler法的计算代码如下，这里由于考虑到对比，针对同一个ODE进行计算，只不过需要利用两种方法，因此可以通过构建一个ODE类，而两个方法可以绑定为实例的方法，并且提供一个可视化方法。

import numpy as np
import matplotlib.pyplot as plt

class Euler_Method_Analysis(object):
    def __init__(self, interval, y0, n, df):
        self.h = (interval[-1] - interval[0]) / n
        self.t = np.arange(interval[0], interval[-1]+self.h, self.h)
        self.y0 = y0
        self.df = df
        self.classical_euler()
        self.modified_euler()
    
    def classical_euler(self):
        w = np.zeros(self.t.shape[0])
        w[0] = self.y0
        for i in np.arange(self.t.shape[0]-1):
            w[i+1] = w[i] + self.h * self.df(self.t[i], w[i])
        self.w_classical = w
        return self.w_classical
    
    def modified_euler(self):
        w = np.zeros(self.t.shape[0])
        w[0] = self.y0
        for i in np.arange(self.t.shape[0] - 1):
            w[i + 1] = w[i] + 0.5 * self.h * (self.df(self.t[i], w[i]) + 
             self.df(self.t[i] + 
                     self.h, w[i] + self.h * self.df(self.t[i], w[i])))
        self.w_modified = w
        return self.w_modified
    
    def view_result(self):
        plt.style.use('ggplot')
        fig, ax = plt.subplots()
        ax.plot(self.t, self.w_classical, lw=2, marker='o')
        ax.plot(self.t, self.w_modified, lw=2, marker='s')
        ax.set_title('Result of $ df = ty + t^3$')
        ax.set_xlabel('t')
        ax.set_ylabel('y')
        label = ['Classical Euler Method', 'Modified Euler Method']
        ax.legend(label, loc='best')
        fig.savefig('resfig.png', dpi=600)
        plt.show()
        
        
# 一下为调用上述类，实现计算
def func(t, y):
    return t * y + t**3

res = Euler_Method_Analysis([0, 1], 1, 10, func)
res.w_classical
res.w_modified
res.view_result()

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！