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

本文主要是介绍东南大学研究生-数值分析上机题（2023）Python 6 常微分方程数值解法，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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

6.1 题目

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

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预测校正）的精度更高；
• 要根据不同的问题选择合适的数值解法，公式的精度越高不代表实际的求解精度越高；
• 常微分方程的数值解法是广泛应用的方法，在以后的工程实践与科研之中会有更多的应用.

这篇关于东南大学研究生-数值分析上机题（2023）Python 6 常微分方程数值解法的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

---