隐式

按方法整理/常微分方程-阿达姆斯隐式方法.py

@@ -0,0 +1,61 @@
+import math
+
+
+def AdamusImplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
+    betas = [
+        [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):
+        print("k超出范围")
+        return None
+    if len(first_ys) != k:
+        print("first_ys(前几个y的值)与k长度不匹配")
+        return None
+    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k)]
+    for i in range(k):
+        result.append((x0 + i * h, first_ys[i]))
+    x1 = x0 + h*k
+    y0 = first_ys[k-1]
+    while x1 < xk:
+        y1_t = y0
+        y1 = y0 + 1
+        for j in range(1000):
+            y1 = y0 + h * sum(betas[k][i+1] * fxys[k-i-1] for i in range(k))+ h * betas[k][0] * fxy(x1, y1_t)
+            if abs(y1 - y1_t) < 1e-14:
+                break
+            y1_t = y1
+
+        # y1 = y0 + delta_y
+        print(f"x={x1}, y={y1}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
+        result.append((x1, y1))
+        fxys.append(fxy(x1, y1))
+        fxys.pop(0)
+        y0 = y1
+        x1 += h
+    return result
+
+
+
+if __name__ == "__main__":
+    # 定义初始值和参数#####################################################################################
+    k = 3                    # 阿达姆斯k步隐式方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P255
+    first_ys = [1, 0.904837418036, 0.818730753078]  # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
+    x0 = 0.0                 # 初始x值
+    y0 = 1.0                 # 初始y值
+    h = 0.1                  # 步长
+    xk = 1.0                 # 最终x值
+    fxy = lambda x, y: -y                   # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
+    fx_real = lambda x: math.exp(-x)        # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
+    # 调用阿达姆斯隐式方法
+    result = AdamusImplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
+    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
+    if result:
+        for xy in result:
+            print(f"(x, y): {xy}")