编写微分方程 dy/dx＝xy, 当 x＝0 时 y＝1+1/5, x 属于 0～3 之间，编写积分程序，包括欧 拉数值积分程序，预报校正数字积分程序、4 阶龙格库塔积分程序，它们的积分步长分别取0.01，0.1, 0.5， 绘制积分结果曲线

根据微分方程 dy/dx=xy，可以进行分离变量得到：

dy/y = x dx

两边同时积分，得到：

ln|y| = x^2/2 + C

其中 C 为常数，代入 x=0，y=1+1/5 可以求得 C = ln(6/5)

因此，原微分方程的通解为：

y = Ce^(x^2/2)，其中 C = 6/5

下面分别编写三种积分程序，并绘制积分结果曲线。

1. 欧拉数值积分程序：

def euler_integration(f, x0, y0, h, n):
x = [x0]
y = [y0]
for i in range(n):
y.append(y[-1] + h * f(x[-1], y[-1]))
x.append(x[-1] + h)
return x, y

def f(x, y):
return x * y

x0 = 0
y0 = 6/5
h1 = 0.01
h2 = 0.1
h3 = 0.5
n = int((3 - x0) / h1)

x1, y1 = euler_integration(f, x0, y0, h1, n)
x2, y2 = euler_integration(f, x0, y0, h2, n)
x3, y3 = euler_integration(f, x0, y0, h3, n)

import matplotlib.pyplot as plt

plt.plot(x1, y1, label='h=0.01')
plt.plot(x2, y2, label='h=0.1')
plt.plot(x3, y3, label='h=0.5')
plt.legend()
plt.show()

2. 预报校正数字积分程序：

def predictor_corrector_integration(f, x0, y0, h, n):
x = [x0]
y = [y0]
for i in range(n):
y_predict = y[-1] + h * f(x[-1], y[-1])
y_correct = y[-1] + h * (f(x[-1], y[-1]) + f(x[-1] + h, y_predict)) / 2
y.append(y_correct)
x.append(x[-1] + h)
return x, y

x0 = 0
y0 = 6/5
h1 = 0.01
h2 = 0.1
h3 = 0.5
n = int((3 - x0) / h1)

x1, y1 = predictor_corrector_integration(f, x0, y0, h1, n)
x2, y2 = predictor_corrector_integration(f, x0, y0, h2, n)
x3, y3 = predictor_corrector_integration(f, x0, y0, h3, n)

plt.plot(x1, y1, label='h=0.01')
plt.plot(x2, y2, label='h=0.1')
plt.plot(x3, y3, label='h=0.5')
plt.legend()
plt.show()

3. 4阶龙格库塔积分程序：

def runge_kutta_integration(f, x0, y0, h, n):
x = [x0]
y = [y0]
for i in range(n):
k1 = h * f(x[-1], y[-1])
k2 = h * f(x[-1] + h/2, y[-1] + k1/2)
k3 = h * f(x[-1] + h/2, y[-1] + k2/2)
k4 = h * f(x[-1] + h, y[-1] + k3)
y.append(y[-1] + (k1 + 2*k2 + 2*k3 + k4) / 6)
x.append(x[-1] + h)
return x, y

x0 = 0
y0 = 6/5
h1 = 0.01
h2 = 0.1
h3 = 0.5
n = int((3 - x0) / h1)

x1, y1 = runge_kutta_integration(f, x0, y0, h1, n)
x2, y2 = runge_kutta_integration(f, x0, y0, h2, n)
x3, y3 = runge_kutta_integration(f, x0, y0, h3, n)

plt.plot(x1, y1, label='h=0.01')
plt.plot(x2, y2, label='h=0.1')
plt.plot(x3, y3, label='h=0.5')
plt.legend()
plt.show()