东南大学研究生-数值分析上机题（2023）Python 6 常微分方程数值解法

最新推荐文章于 2024-09-15 22:31:42 发布

天空的蓝耀

最新推荐文章于 2024-09-15 22:31:42 发布

阅读量636

点赞数 9

文章标签： python

本文链接：https://blog.csdn.net/qq_54716613/article/details/141815203

版权

常微分方程初值问题数值解

6.1 题目

编制RK4方法的通用程序；
编制AB4方法的通用程序（由RK4提供初值）；
编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；
编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；
对于初值问题
$\begin{cases} y'=-x^{2}y^{2}, & 0\leq x \leq 1.5,\\ y(0)=3 & \\ \end{cases}$
取步长 $h = 0.1$ ，应用(1)-(4)中的四种方法进行计算，并将计算结果和精确解 $y(x)=3/(1+x^3)$ 作比较；
通过本上机题，你能得到哪些结论？

6.2 Python源程序

# 定义一阶微分方程  
def y_fxy(x, y):  
    return - (x ** 2) * (y ** 2)  
  
# 定义一阶微分方程的精确解 函数  
def y(x):  
    return 3 / (1 + x ** 3)  
  
# RK4  
def rk4(y0, h_, x0, xi):  
    N = int((xi - x0) / h_)  # 运算的次数  
    x = [x0]  
    y_pdt = [y0]  # 近似解  
    y_real = [y0]  
    err_list = [y0-y0]  
  
    for i in range(N):  
        k1 = y_fxy(x[-1], y_pdt[-1])  
        k2 = y_fxy(x[-1] + 1 / 2 * h_, y_pdt[-1] + 1 / 2 * h_ * k1)  
        k3 = y_fxy(x[-1] + 1 / 2 * h_, y_pdt[-1] + 1 / 2 * h_ * k2)  
        k4 = y_fxy(x[-1] + h_, y_pdt[-1] + h_ * k3)  
        y_pdt.append(y_pdt[-1] + h_ / 6 * (k1 + 2 * k2 + 2 * k3 + k4))  
        x.append(x[-1]+h_)  
        y_real.append(y(x[-1]))  
        err_list.append(y_real[-1] - y_pdt[-1])  
  
    return x, y_pdt, y_real, err_list  
  
# AB4  
def ab4(y0, h_, x0, xi):  
    N = int((xi - x0) / h_)  # 运算的次数  
  
    x, y_pdt, y_real, err_list = rk4(y0, h_, x0, x0 + 3 * h_) # y0给定 y1,y2,y3由RK4得出  
  
    for i in range(3, N):  
        y_pdt.append(y_pdt[-1] + h_ / 24 * \  
                     (55 * y_fxy(x[-1], y_pdt[-1]) - 59 * y_fxy(x[-2], y_pdt[-2]) + \  
                      37 * y_fxy(x[-3], y_pdt[-3]) - 9 * y_fxy(x[-4], y_pdt[-4])))  
        x.append(x[-1]+h_)  
        y_real.append(y(x[-1]))  
        err_list.append(y_real[-1] - y_pdt[-1])  
  
    return x, y_pdt, y_real, err_list  
  
# AB4_AM4预测算法  
def ab4_am4(y0, h_, x0, xi):  
    N = int((xi - x0) / h_)  # 运算的次数  
  
    x, y_pdt, y_real, err_list = rk4(y0, h_, x0, x0 + 3 * h_) # y0给定 y1,y2,y3由RK4得出  
  
    for i in range(3, N):  
        y_pdt.append(y_pdt[-1] + h_ / 24 * \  
                     (55 * y_fxy(x[-1], y_pdt[-1]) - 59 * y_fxy(x[-2], y_pdt[-2]) + \  
                      37 * y_fxy(x[-3], y_pdt[-3]) - 9 * y_fxy(x[-4], y_pdt[-4])))  
        x.append(x[-1] + h_)  
        y_pdt[-1] = y_pdt[-2] + h_ / 24 * \  
                     (9 * y_fxy(x[-1], y_pdt[-1]) + 19 * y_fxy(x[-2], y_pdt[-2]) - \  
                      5 * y_fxy(x[-3], y_pdt[-3]) + y_fxy(x[-4], y_pdt[-4]))  
        y_real.append(y(x[-1]))  
        err_list.append(y_real[-1] - y_pdt[-1])  
  
    return x, y_pdt, y_real, err_list  
  
# 改进的AB4_AM4预测算法  
def plus_ab4_am4(y0, h_, x0, xi):  
    N = int((xi - x0) / h_)  # 运算的次数  
  
    x, y_pdt, y_real, err_list = rk4(y0, h_, x0, x0 + 3 * h_) # y0给定 y1,y2,y3由RK4得出  
  
    for i in range(3, N):  
        y_pdt.append(y_pdt[-1] + h_ / 24 * \  
                     (55 * y_fxy(x[-1], y_pdt[-1]) - 59 * y_fxy(x[-2], y_pdt[-2]) + \  
                      37 * y_fxy(x[-3], y_pdt[-3]) - 9 * y_fxy(x[-4], y_pdt[-4])))  
        x.append(x[-1] + h_)  
        y_c = y_pdt[-2] + h_ / 24 * \  
                     (9 * y_fxy(x[-1], y_pdt[-1]) + 19 * y_fxy(x[-2], y_pdt[-2]) - \  
                      5 * y_fxy(x[-3], y_pdt[-3]) + y_fxy(x[-4], y_pdt[-4]))  
        y_pdt[-1] = 251 / 270 * y_c + 19 / 270 * y_pdt[-1]  
        y_real.append(y(x[-1]))  
        err_list.append(y_real[-1] - y_pdt[-1])  
  
    return x, y_pdt, y_real, err_list  
  
  
def display(x, y_pdt, y_real, err_list, h_, x0, xi):  
    N = int((xi - x0) / h_)  # 运算的次数  
    print("i  xi      yi        y(xi)    y(xi)-yi")  
    for i in range(N):  
        print("{:d} {:.2f} {:.8f} {:.8f} {:.8f}".format\  
                  (i+1, x[i+1], y_pdt[i+1], y_real[i+1], err_list[i+1]))  
  
  
if __name__ == '__main__':  
    y_0 = 3  # 初值  
    h = 0.1  # 步长  
    x_0 = 0  # 区间左端点  
    x_i = 1.5  # 区间右端点  
  
    X, Y_pdt, Y_real, Error = rk4(y_0, h, x_0, x_i)  
    print("RK4:")  
    display(X, Y_pdt, Y_real, Error, h, x_0, x_i)  
    print("RK4整体截断误差：{:.8f}".format(max(list(map(abs, Error)))))  
  
    X, Y_pdt, Y_real, Error = ab4(y_0, h, x_0, x_i)  
    print("AB4:")  
    display(X, Y_pdt, Y_real, Error, h, x_0, x_i)  
    print("AB4整体截断误差：{:.8f}".format(max(list(map(abs, Error)))))  
  
    X, Y_pdt, Y_real, Error = ab4_am4(y_0, h, x_0, x_i)  
    print("AB4-AM4预测校正:")  
    display(X, Y_pdt, Y_real, Error, h, x_0, x_i)  
    print("AB4-AM4预测矫正整体截断误差：{:.8f}".format(max(list(map(abs, Error)))))  
  
    X, Y_pdt, Y_real, Error = plus_ab4_am4(y_0, h, x_0, x_i)  
    print("改进的AB4-AM4预测校正:")  
    display(X, Y_pdt, Y_real, Error, h, x_0, x_i)  
    print("改进的AB4-AM4预测矫正整体截断误差：{:.8f}".format(max(list(map(abs, Error)))))

6.3 程序运行结果

RK4:

RK4:
i  xi      yi        y(xi)    y(xi)-yi
1 0.10 2.99700281 2.99700300 0.00000019
2 0.20 2.97619008 2.97619048 0.00000039
3 0.30 2.92112875 2.92112950 0.00000076
4 0.40 2.81954726 2.81954887 0.00000161
5 0.50 2.66666349 2.66666667 0.00000318
6 0.60 2.46710026 2.46710526 0.00000501
7 0.70 2.23379914 2.23380491 0.00000577
8 0.80 1.98412285 1.98412698 0.00000413
9 0.90 1.73510711 1.73510700 -0.00000012
10 1.00 1.50000581 1.50000000 -0.00000581
11 1.10 1.28701259 1.28700129 -0.00001131
12 1.20 1.09972217 1.09970674 -0.00001542
13 1.30 0.93839746 0.93837973 -0.00001773
14 1.40 0.80130043 0.80128205 -0.00001838
15 1.50 0.68573209 0.68571429 -0.00001780
RK4整体截断误差：0.00001838

AB4:

AB4:
i  xi      yi        y(xi)    y(xi)-yi
1 0.10 2.99700281 2.99700300 0.00000019
2 0.20 2.97619008 2.97619048 0.00000039
3 0.30 2.92112875 2.92112950 0.00000076
4 0.40 2.81838926 2.81954887 0.00115961
5 0.50 2.66467247 2.66666667 0.00199420
6 0.60 2.46520263 2.46710526 0.00190263
7 0.70 2.23307895 2.23380491 0.00072596
8 0.80 1.98495058 1.98412698 -0.00082359
9 0.90 1.73704329 1.73510700 -0.00193629
10 1.00 1.50219455 1.50000000 -0.00219455
11 1.10 1.28876344 1.28700129 -0.00176216
12 1.20 1.10072420 1.09970674 -0.00101746
13 1.30 0.93871050 0.93837973 -0.00033077
14 1.40 0.80113495 0.80128205 0.00014710
15 1.50 0.68533458 0.68571429 0.00037971
AB4整体截断误差：0.00219455

AB4-AM4预测校正:

AB4-AM4预测校正:
i  xi      yi        y(xi)    y(xi)-yi
1 0.10 2.99700281 2.99700300 0.00000019
2 0.20 2.97619008 2.97619048 0.00000039
3 0.30 2.92112875 2.92112950 0.00000076
4 0.40 2.81967843 2.81954887 -0.00012956
5 0.50 2.66687598 2.66666667 -0.00020932
6 0.60 2.46725176 2.46710526 -0.00014650
7 0.70 2.23373141 2.23380491 0.00007350
8 0.80 1.98378670 1.98412698 0.00034028
9 0.90 1.73460744 1.73510700 0.00049956
10 1.00 1.49951594 1.50000000 0.00048406
11 1.10 1.28665714 1.28700129 0.00034415
12 1.20 1.09953315 1.09970674 0.00017360
13 1.30 0.93834252 0.93837973 0.00003721
14 1.40 0.80132737 0.80128205 -0.00004532
15 1.50 0.68579611 0.68571429 -0.00008183
AB4-AM4预测校正整体截断误差：0.00049956

改进的AB4-AM4预测校正:

改进的AB4-AM4预测校正:
i  xi      yi        y(xi)    y(xi)-yi
1 0.10 2.99700281 2.99700300 0.00000019
2 0.20 2.97619008 2.97619048 0.00000039
3 0.30 2.92112875 2.92112950 0.00000076
4 0.40 2.81958771 2.81954887 -0.00003884
5 0.50 2.66671285 2.66666667 -0.00004619
6 0.60 2.46709703 2.46710526 0.00000823
7 0.70 2.23368249 2.23380491 0.00012242
8 0.80 1.98388468 1.98412698 0.00024230
9 0.90 1.73480801 1.73510700 0.00029899
10 1.00 1.49973191 1.50000000 0.00026809
11 1.10 1.28682068 1.28700129 0.00018061
12 1.20 1.09962178 1.09970674 0.00008496
13 1.30 0.93836732 0.93837973 0.00001242
14 1.40 0.80131135 0.80128205 -0.00002930
15 1.50 0.68576045 0.68571429 -0.00004616
改进的AB4-AM4预测校正整体截断误差：0.00029899

6.4 总结感悟

根据数值分析理论推导的结果，RK4、AB4、AB4-AM4预测校正具有4阶精度，而改进的AB4-AM4预测校正具有5阶精度，但是对于该问题来说，比较四种常微分方程数值解法在 $[0.1, 1.5]$ 上的整体截断误差，则是RK4<改进的AB4-AM4预测校正<AB4-AM4预测校正<AB4，RK4（单步法）的精度要比多步法（AB4、AB4-AM4预测校正、改进的AB4-AM4预测校正）的精度更高；
要根据不同的问题选择合适的数值解法，公式的精度越高不代表实际的求解精度越高；
常微分方程的数值解法是广泛应用的方法，在以后的工程实践与科研之中会有更多的应用.