快捷搜索: 王者荣耀 脱发

欧拉方法python代码实现

如有不对,请批评指正!

1. 显式Euler和隐式Euler

2. 改进的Euler

3. 龙格-库塔(Runge-Kutta)

4. 举个栗子

代码仅为了说明计算过程,并未进行优化!

import math
import matplotlib.pyplot as plt

y_e = []  # 显式
y_i = []  # 隐式
y_pro = []  # 改进
y_p = []  # 改进
y_rk = []  # 龙格-库塔
y_c = []  # 解析解
x = []
h = 0.1  # 步长
num = 20


def fun(x0, y0):
    return -2 * x0 / y0 + y0


def Euler():
    print({:>3s}{:>8s}{:>9s}{:>12s}{:>9s}{:>9s}.format(x, y_e, y_i, y_pro, y_rk, y_c))
    for i in range(num):
        # 显示Euler
        y_e.append(round(y_e[i] + h * fun(x[i], y_e[i]), 5))

        # 隐式Euler
        y_i1 = y_i[i] + h * fun(x[i], y_i[i])
        y_i2 = y_i[i] + h * fun(x[i + 1], y_i1)
        while abs(y_i2 - y_i1) > 1e-6:
            y_i1 = y_i2
            y_i2 = y_i[i] + h * fun(x[i + 1], y_i1)
        y_i.append(y_i2)

        # 改进Euler
        y_p.append(y_pro[i] + h * fun(x[i], y_pro[i]))
        y_pro.append(round(y_pro[i] + h / 2 * (fun(x[i], y_pro[i]) + fun(x[i + 1], y_p[i])), 5))

        # 龙格-库塔
        k1 = fun(x[i], y_rk[i])
        k2 = fun(x[i] + h / 2, y_rk[i] + h / 2 * k1)
        k3 = fun(x[i] + h / 2, y_rk[i] + h / 2 * k2)
        k4 = fun(x[i] + h, y_rk[i] + h * k3)
        y_rk.append(round(y_rk[i] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4), 5))

        # 解析解
        y_c.append(round(math.sqrt(1 + 2 * x[i + 1]), 5))

        print({}{:>10f}{:>10f}{:>10f}{:>10f}{:>10f}
.format(
            round(x[i + 1], 1),
            y_e[i + 1],
            y_i[i + 1],
            y_pro[i + 1],
            y_rk[i + 1],
            y_c[i + 1]))

    plt.plot(x, y_e, label=Euler_Explicit, color=r, marker=+)
    plt.plot(x, y_i, label=Euler_Implicit, color=c, marker=>)
    plt.plot(x, y_pro, label=Euler_Improvement, color=b, marker=3)
    plt.plot(x, y_rk, label=Euler_RungeKutta, color=m, marker=d)
    plt.plot(x, y_c, label=Analytical solution, color=g, marker=o)
    plt.legend()
    plt.show()


if __name__ == __main__:
    # 初始值
    x.append(0.0)
    y_e.append(1)  # 显式
    y_i.append(1)  # 隐式
    y_pro.append(1)  # 改进
    y_rk.append(1)  # 龙格-库塔
    y_c.append(1)  # 解析解
    for i in range(num):
        x.append(round(x[0] + (i + 1) * h, 1))
    Euler()

5. 运行结果

num = 6时

num = 20

参考:

经验分享 程序员 微信小程序 职场和发展