12-11：求微分欧拉法

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

def f(x,y):
    return y-(2*x/y)

def Y(x):
    return np.sqrt(1+2*x)

def Euler(a,b,h,y0):
    x = np.arange(a,b+h,h)
    yy = [y0]
    for i in x:
        y = y0 + h*f(i,y0)
        yy.append(y)
        y0 = y
    return x,np.array(yy[:-1])

def EulerPlus(a,b,h,y0):
    x = np.arange(a,b+h,h)
    yy = [y0]
    for i in range(len(x)-1):
        y = y0 + h*f(x[i],y0)

        ynext = y0 + (h/2)*(f(x[i],y0)+f(x[i+1],y))

        yy.append(ynext)
        y0 = ynext
    re = np.array(yy)
    return x,re

a,b,h,y0 = 0,3,0.1,1
x,yk = Euler(a,b,h,y0)

plt.scatter(x,yk,c='r')
plt.plot(x,Y(x))

output_2_1.png

改进欧拉法

a,b,h,y0 = 0,3,0.1,1
x,re = EulerPlus(a,b,h,y0)

plt.figure(figsize=(8,6))
plt.plot(x,Y(x))
plt.scatter(x,yk)
plt.scatter(x,re,c='r')
plt.legend(['真实值','显示欧拉法','改进欧拉法'])
plt.grid()

output_4_0.png