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
