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

2024-09-02 15:44
文章标签 python 分析 2023 研究生 上机 数值 微分方程 解法 东南大学

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

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

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

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


http://www.chinasem.cn/article/1130375。

相关文章

python: 多模块(.py)中全局变量的导入

python: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别： 当你在一个模块（文件）中使用 global 声明变量时，这个变量只在该模块的全局命名空
阅读更多...
性能分析之MySQL索引实战案例

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题，本文将一起从代码层去分析为什么没有建立索引？ 开源ERP项目地址：https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置
阅读更多...
【Python编程】Linux创建虚拟环境并配置与notebook相连接

【Python编程】Linux创建虚拟环境并配置与notebook相连接

1.创建 使用 venv 创建虚拟环境。例如，在当前目录下创建一个名为 myenv 的虚拟环境： python3 -m venv myenv 2.激活 激活虚拟环境使其成为当前终端会话的活动环境。运行： source myenv/bin/activate 3.与notebook连接 在虚拟环境中，使用 pip 安装 Jupyter 和 ipykernel： pip instal
阅读更多...
poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大，用原来的O（n^2）会爆。 小白上的O(n^2)解法代码：TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO
阅读更多...
【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

引言 高斯过程（Gaussian Process，简称GP）是一种概率模型，用于描述一组随机变量的联合概率分布，其中任何一个有限维度的子集都具有高斯分布 文章目录 引言一、高斯过程1.1 基本定义1.1.1 随机过程1.1.2 高斯分布 1.2 高斯过程的特性1.2.1 联合高斯性1.2.2 均值函数1.2.3 协方差函数（或核函数） 1.3 核函数1.4 高斯过程回归（Gauss
阅读更多...
【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络（1）sklearn

系列文章目录 监督学习：参数方法 【学习笔记】 陈强-机器学习-Python-Ch4 线性回归 【学习笔记】 陈强-机器学习-Python-Ch5 逻辑回归 【课后题练习】 陈强-机器学习-Python-Ch5 逻辑回归（SAheart.csv） 【学习笔记】 陈强-机器学习-Python-Ch6 多项逻辑回归 【学习笔记 及 课后题练习】 陈强-机器学习-Python-Ch7 判别分析 【学
阅读更多...
nudepy，一个有趣的 Python 库！

nudepy，一个有趣的 Python 库！

更多资料获取 📚 个人网站：ipengtao.com 大家好，今天为大家分享一个有趣的 Python 库 - nudepy。 Github地址：https://github.com/hhatto/nude.py 在图像处理和计算机视觉应用中，检测图像中的不适当内容（例如裸露图像）是一个重要的任务。nudepy 是一个基于 Python 的库，专门用于检测图像中的不适当内容。该
阅读更多...
SWAP作物生长模型安装教程、数据制备、敏感性分析、气候变化影响、R模型敏感性分析与贝叶斯优化、Fortran源代码分析、气候数据降尺度与变化影响分析

SWAP作物生长模型安装教程、数据制备、敏感性分析、气候变化影响、R模型敏感性分析与贝叶斯优化、Fortran源代码分析、气候数据降尺度与变化影响分析

查看原文>>>全流程SWAP农业模型数据制备、敏感性分析及气候变化影响实践技术应用 SWAP模型是由荷兰瓦赫宁根大学开发的先进农作物模型，它综合考虑了土壤-水分-大气以及植被间的相互作用；是一种描述作物生长过程的一种机理性作物生长模型。它不但运用Richard方程，使其能够精确的模拟土壤中水分的运动，而且耦合了WOFOST作物模型使作物的生长描述更为科学。 本文让更多的科研人员和农业工作者
阅读更多...
MOLE 2.5 分析分子通道和孔隙

MOLE 2.5 分析分子通道和孔隙

软件介绍 生物大分子通道和孔隙在生物学中发挥着重要作用，例如在分子识别和酶底物特异性方面。 我们介绍了一种名为 MOLE 2.5 的高级软件工具，该工具旨在分析分子通道和孔隙。 与其他可用软件工具的基准测试表明，MOLE 2.5 相比更快、更强大、功能更丰富。作为一项新功能，MOLE 2.5 可以估算已识别通道的物理化学性质。 软件下载 https://pan.quark.cn/s/57
阅读更多...
pip-tools：打造可重复、可控的 Python 开发环境，解决依赖关系，让代码更稳定

pip-tools：打造可重复、可控的 Python 开发环境，解决依赖关系，让代码更稳定

在 Python 开发中，管理依赖关系是一项繁琐且容易出错的任务。手动更新依赖版本、处理冲突、确保一致性等等，都可能让开发者感到头疼。而 pip-tools 为开发者提供了一套稳定可靠的解决方案。 什么是 pip-tools？ pip-tools 是一组命令行工具，旨在简化 Python 依赖关系的管理，确保项目环境的稳定性和可重复性。它主要包含两个核心工具：pip-compile 和 pip
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2