当前位置： 首页 > article >正文

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

article 2025/1/22 8:09:31

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

6.1 题目

  1. 编制RK4方法的通用程序；
  2. 编制AB4方法的通用程序（由RK4提供初值）；
  3. 编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；
  4. 编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；
  5. 对于初值问题
    { y ′ = − x 2 y 2 , 0 ≤ x ≤ 1.5 , y ( 0 ) = 3 \begin{cases} y'=-x^{2}y^{2}, & 0\leq x \leq 1.5,\\ y(0)=3 & \\ \end{cases} {y′=−x2y2,y(0)=3​0≤x≤1.5,​
    取步长 h = 0.1 h=0.1 h=0.1，应用(1)-(4)中的四种方法进行计算，并将计算结果和精确解 y ( x ) = 3 / ( 1 + x 3 ) y(x)=3/(1+x^3) y(x)=3/(1+x3)作比较；
  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 ] [0.1,1.5] [0.1,1.5]上的整体截断误差，则是RK4<改进的AB4-AM4预测校正<AB4-AM4预测校正<AB4，RK4（单步法）的精度要比多步法（AB4、AB4-AM4预测校正、改进的AB4-AM4预测校正）的精度更高；
  • 要根据不同的问题选择合适的数值解法，公式的精度越高不代表实际的求解精度越高；
  • 常微分方程的数值解法是广泛应用的方法，在以后的工程实践与科研之中会有更多的应用.

http://www.kler.cn/a/287118.html

相关文章：

  • 玉米植物结构受乙烯生物合成基因 ZmACS7 的调控
  • 2025CSP-J 冲刺训练(3)：前缀和差分
  • elasticsearch 数据导出/导入
  • Windows FileZila Server共享电脑文件夹 映射21端口外网连接
  • .Net Core微服务入门全纪录（五）——Ocelot-API网关（下）
  • Docker使用 使用Dockerfile来创建镜像
  • CSS3中display显示属性
  • ChatGPT 3.5/4.0新手使用手册
  • SQL语句中模糊匹配LIKE和正则表达式之间有什么相同点和不同点
  • compser好用镜像
  • 基于Python的Flask框架实战全流程从新建到部署【2】
  • Atlas阿特拉斯wordpress主题
  • 编译与运行
  • C++：二叉树进阶
  • Ubuntu18.04 下安装CUDA
  • 企业IT服务管理（ITSM）的实践与探索
  • 永磁同步电机高性能控制算法（15）——自抗扰控制电流环的详细解读/最新TIE论文解析
  • uni-app 获取当前位置的经纬度以及地址信息
  • NRK3301语音识别芯片免联网实现智能空气炸锅
  • 【系统架构设计师】论文：论混合软件架构的设计
  • NLP从零开始------16.文本中阶处理之序列到序列模型（1）
  • sql-libs第三关详细解答
  • 1412312
  • shell编程--数组
  • Leetcode 77. 组合 组合型回溯 C++实现
  • 【STM32】红外遥控