# 常微分方程初值问题的数值方法 Python实现

#### 原理：

（1）利用Euler方法和改进的Euler方法求解初值问题
（2）利用Runge-Kutta方法求解初值问题

import math

def f1(x, y):
if x == 0:
return 0
else:
return ((4 * x) / y) - (x * y)
def f2(x, y):
return x ** 2 - y ** 2

def Euler(a, b, alpha, N):
h = (b - a) / N
x = a
y = alpha
out = [(x, y, y, y-y)]
print ("%-10s %-10s %-10s %-10s"%("Xn","Yn","Y(Xn)","Y(Xn)-Yn"))
print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,y,y-y))
for n in range(1, N + 1):
y = y + h * f1(x,y)
x = x + h
Y = (4 + 5 * (math.e ** -(x ** 2))) ** 0.5
out.append((x, y, Y, y-Y))
print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,Y,y-Y))
return out

def EulerPlus(a, b, alpha, N):
h = (b - a) / N
x = a
y = alpha
print ("%-10s %-10s %-10s %-10s"%("Xn","Yn","Y(Xn)","Y(Xn)-Yn"))
print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,y,y-y))
out = [(x, y, y, y-y)]
for n in range(1, N + 1):
K1 = f1(x, y)
K2 = f1((x + h), (y + h * K1))
y = y + h * (K1 + K2) / 2
x = x + h
Y = (4 + 5 * (math.e ** -(x ** 2))) ** 0.5
out.append((x, y, Y, y-Y))
print ("%-10.6f %-10.6f %-10.6f %-10.6f"%(x,y,Y,y-Y))
return out

def RungeKutta(a, b, alpha, N):
h = (b - a) / N
x = a
y = alpha
out = [(x,y)]
print ("%-10s %-10s"%("Xn","Yn"))
print ("%-10.6f %-10.6f"%(x,y))
for n in range(1, N + 1):
K1 = f2(x, y)
K2 = f2((x + h / 2),(y + (h * K1) / 2))
K3 = f2((x + h / 2),(y + (h * K2) / 2))
K4 = f2((x + h),(y + h * K3))
y = y + h * (K1 + 2 * K2 + 2 * K3 + K4) / 6
x = x + h
out.append((x,y))
print ("%-10.6f %-10.6f"%(x,y))
return out

Euler(0,2,3,20)
Euler(0,2,3,10)
Euler(0,2,3,5)
EulerPlus(0,2,3,20)
EulerPlus(0,2,3,10)
EulerPlus(0,2,3,5)
RungeKutta(-1,0,0,10)
RungeKutta(-1,0,0,5)


02-02 2万+

04-21 8132
07-06 3654
09-24 473
09-11
07-17 568
04-27 434
09-11 4355