65 lines
2.6 KiB
Python
65 lines
2.6 KiB
Python
import math
|
||
|
||
|
||
def AdamusFix(k,first_ys,x0,y0,h,xk,fxy,fx_real):
|
||
betas_1 = [
|
||
[1],
|
||
[3/2, -1/2],
|
||
[23/12, -16/12, 5/12],
|
||
[55/24, -59/24, 37/24, -9/24],
|
||
[1901/720, -2774/720, 2626/720, -1274/720, 251/720],
|
||
[4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
|
||
]
|
||
# Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
|
||
betas_2 = [
|
||
[1],
|
||
[1/2, 1/2],
|
||
[5/12, 8/12, -1/12],
|
||
[9/24, 19/24, -5/24, 1/24],
|
||
[251/720, 646/720, -261/720, 106/720, -19/720],
|
||
[475/1440, 1427/1440, -798/1440, 482/1440, -173/1440, 27/1440]
|
||
]
|
||
result = []
|
||
if k<0 or k >= len(betas_1):
|
||
print("k超出范围")
|
||
return None
|
||
if len(first_ys) != k + 1:
|
||
print("first_ys(前几个y的值)与k长度不匹配")
|
||
return None
|
||
fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
|
||
for i in range(k + 1):
|
||
result.append((x0 + i * h, first_ys[i]))
|
||
x1 = x0 + h*k
|
||
y0 = first_ys[k]
|
||
while x1 < xk:
|
||
x1 += h
|
||
y1_guess = y0 + h * sum(betas_1[k][i] * fxys[k-i] for i in range(k + 1))
|
||
fxy_guess = fxy(x1, y1_guess)
|
||
fxys.append(fxy_guess)
|
||
fxys.pop(0)
|
||
delta_y = h * sum(betas_2[k][i] * fxys[k-i] for i in range(k + 1))
|
||
y1 = y0 + delta_y
|
||
print(f"x={x1}, y={y1}, y_guess={y1_guess}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
|
||
result.append((x1, y1))
|
||
fxys[k] = fxy(x1, y1)
|
||
y0 = y1
|
||
return result
|
||
|
||
|
||
|
||
if __name__ == "__main__":
|
||
# 定义初始值和参数#####################################################################################
|
||
k = 3 # 阿达姆斯k+1步校正方法，精度为k+1阶，最常用k=3，其他阶数我没试过 P253
|
||
first_ys = [1, 1.0954, 1.1832, 1.2649] # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
|
||
x0 = 0.0 # 初始x值
|
||
y0 = 1.0 # 初始y值
|
||
h = 0.1 # 步长
|
||
xk = 1.0 # 最终x值
|
||
fxy = lambda x, y: y - 2*x/y # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
|
||
fx_real = lambda x: 0 # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
|
||
# 调用阿达姆斯显式方法
|
||
result = AdamusFix(k, first_ys, x0, y0, h, xk, fxy, fx_real)
|
||
print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
|
||
if result:
|
||
for xy in result:
|
||
print(f"(x, y): {xy}") |