高教杯数学建模A题程序设计要点与思路

news/2025/2/16 6:20:51/
  • 2023 年是我最后一次参加 高教杯大学生数学建模竞赛 以后不会再参加了（大四参加意义不太，研究生有研究生的数学建模大赛）
    • 很遗憾 由于各种原因 我们没有能够完成赛题
    • 2022 年 美赛 2022年 Mathor Cup 2022 年国赛 2022 亚太杯 2023年 美赛 2023年 国赛
    • 我和我的朋友一共参加了6次比赛
    • 6次比赛 我交到了很好的朋友 
    • 然鹅 成绩比较惨淡 S*1 H*1 省三一次
    • 唉，数模啊数模，很难说出口啊
    • 今年本来特别有信心的，但是，唉，。。。A题，我的痛
  • 为了研究A题
    • 我们先后学习了
      • 数理方程
      • 数值分析
      • 数值仿真
      • 优化算法
      • 稳定性的理论
      • ......
    • 就是想在2023 年A题大展身手，可惜啊可惜
      • 光学 我的痛
      • 真的很难过 心里一直算算的 感觉对不起朋友
      • 确实是我自己不行
      • 自己不行啊
    • 暑假刚做过手术
    • 暑假手术之前去了一些大学
      • 一个优营也没有拿到，毕竟年龄有点小，基础不扎实
    • 本来计划三年物理系毕业然后转强基计划然后深造的，因为我真的修完了课程。
      • 但是计划有变啊。三年毕竟不够扎实
      • 现在去选了 核与粒子物理 欢迎交流
    • 最后一次参加数学建模是2023UPC 好像是11月第一个周末
      • 加油 
  • 本文主要介绍这几个极其重要的程序设计思路，欢迎大家参考，引用
    • 真的，毕竟百度也是可以引用的嘛
    • 不引用我绝对不会在意的

数值常微分

参考书目

  • 微分方程的数值解法与程序实现  华冬英 李祥贵
  • 量子物理学中的常用算法与程序 井孝功 赵永芳 蒿凤有 
  • 稳定性的理论，方法和应用
  • 微分方程 动力系统与混沌导论

简介

  • 很多同学学习数值分析，都是 高斯消元 牛顿迭代 然后...
    • 譬如我写过的一篇博客

Python牛顿迭代法的应用

追赶法(Thomas) 雅克比迭代(Jacobi) 高斯迭代(Gauss) 的C++实现

  • 其实我不建议这样，为什么呢？
    • 函数零点（非线性方程组的解）的求法非常重要，是现代科学计算的基础性内容，但是，这并不意味着我们必须从这些知识学起，原因有2
      • 1.简单的方法收敛条件严格，具体情况需要具体分析，大概率需要更强的算法
      • 2.已经有封装好的先进方法，不应该在基础方法上浪费时间
        • 包括程序设计的时间
        • 包括程序运行时，因为方法优化不得当导致的计算增加时间
  • 数值常微分方程求解包括下面两类
    • 边值问题的求解（本征值问题的求解）
    • 初值问题的求解
      • 辛算法

边值问题

  • 这个问题非常常见
  • 学习过数学物理方程的同学或者泛函分析的课程对此一定印象深刻
  • 我写过的一些博客有

变分原理与边值问题的计算机处理

计算物理专题：双向打靶法解决本征值问题

计算物理专题：有限差分法解决本征值问题

Numerov算法解一维无限深势阱的问题 (含量子力学导论)

  • 你可以使用matlab Pdetool 辅助求解 这也是非常棒的选择
  • 具体的代码就不放在这里了，因为这个难度比较低，关键是验证的难度很低，图一画就知道算得对不对了，也是暴力求解可以解决的问题，完全没有压力

初值问题

  • 这个问题是最常见的问题 譬如2022 年高教杯A题
  • 如果读者对 动力系统有一定了解的话 这个问题可以做的非常漂亮
  • 我个人认为应当从动力系统学起，这样的话分析解的稳定性就会更完整更合理，而不是放一个试验方程的小招就结束了
  • 初值问题的解法分为两类
    • 隐式解法
    • 显式解法
  • 简单得说，显式解法指的是解的当前步完全由解的前几步决定
  • 而隐式解法指的是
    • 只能得到以当前解为未知数的一个方程，必须求解这个方程才能得到当前解
  • 这也就是为什么我不建议花太多时间在高斯迭代，牛顿迭代上的原因。
    • 对于单个的函数，分析算法的稳定是方便的，但是实际问题中我们更关注的是大型方程组的求解，他们往往是非线性的，耦合的，难以分析的，这时解法的稳定性分析难度会非常得大，不如采取：软件+验证的方法
  • 我写过的相关的一些博客我列在这里，具体的代码我就不反复引用了

Python数值分析案例01--------四阶龙格库塔法解抛体运动

常微分方程的龙格库塔显式与隐式解法

  • 尤其是第二篇文章，非常值得好好阅读，
    • 直接引用相关的方法即可，均经过了实践的检验
    • 包含了一个并行加速的例子，可能需要先学习一下multiprocessing 
  • 提示：每次求解必须有稳定性的分析，没有稳定性的分析，你的解的意义何在？可信度何在？

哈密顿系统的辛算法

  • 这是很出彩的一个地方
  • 如果有同学能在数学建模竞赛中实现这一点，我想是极其棒的一件事情，而且我从未见过有很好的论文在本科生竞赛中实现这一点，这是很不现代的
  • 在这里我就不赘述他的优点在什么地方了，较为复杂，公式也很多，
    • 但是只说一句话吧，这个方法算得出彩，就行了
  • 这是我写的一篇博客，可以看看，这非常的基础，你需要很好阅读参考书才能很好得应用他

自洽可分的哈密顿系统的辛算法

参考文献

  • 量子系统的辛算法 丁培柱

数值积分

  • 数值积分常常在数学建模竞赛的某一处出现，还是非常重要的

蒙特卡洛积分方法

  • 蒙特卡洛积分处理高维积分的效果会更优秀一点，在我很近的一篇文章中提到了对比，但只是很粗略的一笔蒙特卡洛方法的数学基础-1
  • 所以如果只是处于装一下的需要的话，完全没有必要写上去

近似积分方法

  • 大家常用的 矩形近似法 梯形近似法 辛普森近似 都属于这一类
  • 我想，如果是已知了具体函数形式的话，使用外推法会更好一些

计算物理专题：高维Romberg数值积分方法

  • 但很多时候，我们只有每个点的函数值，那么还是使用辛普森方法或者更高阶的方法要好

中心差分式

  • 中心差分式是大家都很熟悉的东西了
  • 我需要提到的一点是，你选取的方法的精度尽量要高一些，譬如
    • 目标是 二阶近似
    • 某一步需要用到中心差分值的方法是四阶近似的
    • 那么你的中心差分式精度的选择应该不低于四阶才好

数值偏微分

  • 这也是常考的问题
  • 但是数值偏微分方程的求解极其得复杂，稳定性特别难以分析，我举一个 用迎风法求解抛物型方程的例子
import numpy as np
import matplotlib.pyplot as plt#\frac{\partial u}{\partial t} = \lambda(
#\frac{\partial^2 u}{\partial x^2} +
#\frac{\partial^2 u}{\partial y^2}
#) + q_vclass projequation2D():def __init__(self, Lambda=1, qv = lambda t, x, y:0,X_start=0, X_end=1,Y_start=0, Y_end=1,T_start=0, T_end=1,dx = 0.05,dy = 0.05,dt = 0.01):#c = lambda x,y,t:(?)self.Lambda = Lambdaself.qv = qvself.dx = dxself.dy = dyself.dt = dtself.X_start = X_startself.Y_start = Y_startself.T_start = T_startself.X_end = X_endself.Y_end = Y_endself.T_end = T_endself.X = np.arange(X_start, X_end, dx)self.Y = np.arange(Y_start, Y_end, dy)self.T = np.arange(T_start, T_end, dt)self.results = np.zeros((len(self.T), len(self.X), len(self.Y)))self.startresults()def initcondition_0(self):X_num = len(self.X)Y_num = len(self.Y)test_fun = lambda x,y: np.sin(5*x)*np.cos(5*y)for ix in range(X_num):for iy in range(Y_num):self.results[0][ix][iy] == test_fun(self.X[ix], self.Y[iy])def boundarycondition_left(self):Y_num = len(self.Y)T_num = len(self.T)for it in range(T_num):for iy in range(Y_num):self.results[it][0][iy] = 25def boundarycondition_right(self): Y_num = len(self.Y)T_num = len(self.T)for it in range(T_num):for iy in range(Y_num):self.results[it][-1][iy] = 25def boundarycondition_up(self):X_num = len(self.X)T_num = len(self.T)for it in range(T_num):for iy in range(X_num):self.results[it][iy][-1] = 25def boundarycondition_down(self):X_num = len(self.X)T_num = len(self.T)for it in range(T_num):for iy in range(X_num):self.results[it][iy][0] = 25def startresults(self):self.initcondition_0()self.boundarycondition_left()self.boundarycondition_right()self.boundarycondition_up()self.boundarycondition_down()def Upwind3P(self):for Ti in range(1, len(self.T)):for Xi in range(1, len(self.X)-1):for Yi in range(1, len(self.Y)-1):rx = self.Lambda * self.dt/(self.dx)**2ry = self.Lambda * self.dt/(self.dy)**2self.results[Ti][Xi][Yi] = \rx * (self.results[Ti-1][Xi-1][Yi]+self.results[Ti-1][Xi+1][Yi])+\2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]+\ry * (self.results[Ti-1][Xi][Yi-1]+self.results[Ti-1][Xi][Yi+1])-\self.results[Ti-2][Xi][Yi]return self.resultsdef show_wave(self):fig = plt.figure(figsize=(8,8))ax1 = fig.add_subplot(111, projection='3d')x, y = np.meshgrid(self.X, self.Y)for time in range(len(self.T)):ax1.plot_surface(x, y, self.results[time, :, :],rstride=2, cstride=2, cmap='rainbow')plt.title("time:" + str(self.T[time]))plt.pause(0.5)plt.cla()c = projequation2D()
u = c.Upwind3P()
c.show_wave()
  • 再来一个解波动方程的例子

import numpy as np
import matplotlib.pyplot as plt#\frac{\partial^2 u}{\partial t^2} = c^2 (
#\frac{\partial^2 u}{\partial x^2} +
#\frac{\partial^2 u}{\partial y^2}
#)class waveequation2D():def __init__(self, c,X_start=0, X_end=1,Y_start=0, Y_end=1,T_start=0, T_end=1.01,dx = 0.01,dy = 0.01,dt = 0.01):#c = lambda x,y,t:(?)self.c = cself.dx = dxself.dy = dyself.dt = dtself.X_start = X_startself.Y_start = Y_startself.T_start = T_startself.X_end = X_endself.Y_end = Y_endself.T_end = T_endself.X = np.arange(X_start, X_end, dx)self.Y = np.arange(Y_start, Y_end, dy)self.T = np.arange(T_start, T_end, dt)self.results = np.zeros((len(self.T), len(self.X), len(self.Y)))self.startresults()def initcondition_0(self, x, y):return np.sin(10*x)def initcondition_1(self, x, y):return np.sin(10*x)def boundarycondition_left(self, t, x, y):return (np.sin(10*y))[0]def boundarycondition_right(self, t, x, y): return (np.sin(10*y))[0]def boundarycondition_up(self, t, x, y):return (np.sin(10*x))[0]def boundarycondition_down(self, t, x, y):return (np.sin(10*x))[0]def startresults(self):x, y = np.meshgrid(self.X, self.Y)self.results[0] = self.initcondition_0(x, y)self.results[1] = self.initcondition_1(x, y)t, x, y = np.meshgrid(self.T, self.X, self.Y)self.results[:, 0, :]  = self.boundarycondition_left(t, self.X_start, y)self.results[:, -1, :] = self.boundarycondition_right(t, self.X_end, y)self.results[:, :, -1] = self.boundarycondition_up(t, x, self.Y_end)self.results[:, :, 0]  = self.boundarycondition_down(t, x, self.Y_start)def test_stability(self, x, y, t):test = 4*self.dt**2*self.c(x, y, t)**2/(self.dx**2 + self.dy**2)if test<=1:return Trueelse:print("unstable in x:", x, "y:", y, "t:", t)return Falsedef Upwind3P(self):for Ti in range(2, len(self.T)):for Xi in range(1, len(self.X)-1):for Yi in range(1, len(self.Y)-1):rx = self.c(self.X[Xi], self.Y[Yi], self.T[Ti])**2 * self.dt**2/self.dx**2ry = self.c(self.X[Xi], self.Y[Yi], self.T[Ti])**2 * self.dt**2/self.dy**2self.results[Ti][Xi][Yi] = \rx * (self.results[Ti-1][Xi-1][Yi]+self.results[Ti-1][Xi+1][Yi])+\2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]+\ry * (self.results[Ti-1][Xi][Yi-1]+self.results[Ti-1][Xi][Yi+1])-\self.results[Ti-2][Xi][Yi]
##                    self.test_stability(self.X[Xi],self.Y[Yi],self.T[Ti])                    return self.resultsdef show_wave(self):fig = plt.figure(figsize=(8,12))ax1 = fig.add_subplot(121, projection='3d')ax2 = fig.add_subplot(122, projection='3d')x, y = np.meshgrid(self.X, self.Y)for time in range(len(self.T)):ax1.plot_surface(x, y, self.results[time, :, :],rstride=2, cstride=2, cmap='rainbow')ax2.plot_wireframe(x, y, self.results[time, :, :],rstride=2, cstride=2, linewidth=1, cmap='rainbow')plt.title("time:" + str(self.T[time]))plt.pause(0.01)plt.cla()def test_fun(x, y, t):return 1c = waveequation2D(c=test_fun)
u = c.Upwind3P()
c.show_wave()

  • 所以一般建议使用Matlab 中的pdetool 求解
  • 这是我写的两篇博客，可以参考

典型的偏微分方程数值解法

二维Poisson方程五点差分格式与Python实现

  • 一定要仔细阅读有关书籍，否则很难做出比较好的成果

参考书目

  • 偏微分方程的数值解法（第三版） 陆金甫
  • 特殊函数概论 王竹溪
  • 数学物理方法与仿真（第三版） 杨华军
  • 科学计算中的偏微分方程有限差分法  张文生

优化算法

  • 这是是数学建模的核心，优化求解
  • 一般而言，很多问题可以做凸优化
    • 但是，数学建模中的优化问题往往不能...
    • 所以，要么二分法全局搜索，要么智能求解吧
  • 最近会写一些相关的博客，就不具体演示了。也可以去我的博客下面找......

  • 写完博客，心情舒畅不少


http://www.ppmy.cn/news/1118343.html

相关文章

ELFK之zookeeper+kafka

ELFK之zookeeper+kafka

目录 kafkazookeeper的系统架构 Zookeeper 一、zookeeper概述 二、zookeeper特点 三、zookeeper选举机制 四、应用场景 五、zookeeper实验实例 Kafka 一、概述 为什么需要消息队列(MQ) 使用消息队列的好处 消息队列的两种模式 Kafka 定义 二、Kafka 的特性 三、Ka…
阅读更多...
Springboot整合jdbc和Mybatis

Springboot整合jdbc和Mybatis

目录 整合jdbc 1. 新建项目 2. 编写yaml配置文件连接数据库 3. 测试类 使用原生的jdbcTemplate进行访问测试 使用Druid连接池 1. 添加类型 2. 初始化连接池 3. 编写config类 配置Druid数据源监视 整合Mybatis 1. 导入依赖 2. 编写mapper接口 3. 编写实体类 4. 编…
阅读更多...
如何用ate自动测试设备对DC-DC电源模块负载调整率进行测试？

如何用ate自动测试设备对DC-DC电源模块负载调整率进行测试？

电源模块负载调整率测试是功能测试之一&#xff0c;是电源模块非常重要的一项指标&#xff0c;它的大小直接影响着电源模块的整体质量。因此使用ate自动测试设备对DC-DC电源模块负载调整率进行测试是制造生产过程中至关重要的一环。 电源模块负载调整率计算公式&#xff1a; 负…
阅读更多...
Python基础学习笔记3

Python基础学习笔记3

深度学习实践 深度学习离不开编程 深度学习离不开数学分析&#xff08;高等数学&#xff09;、线性代数、概率论等知识&#xff0c;更离不开以编程为核心的动手实践。 Python编程语言 无论是在机器学习还是深度学习中&#xff0c;Python已经成为主导性的编程语言。而且&…
阅读更多...
面试官：ES6中新增的Set、Map两种数据结构怎么理解?

面试官：ES6中新增的Set、Map两种数据结构怎么理解?

&#x1f3ac; 岸边的风&#xff1a;个人主页 &#x1f525; 个人专栏 :《 VUE 》 《 javaScript 》 ⛺️ 生活的理想&#xff0c;就是为了理想的生活 ! 目录 一、Set 增删改查 add() delete() has() clear() 遍历 二、Map 增删改查 size set() get() has() del…
阅读更多...
gateway之断言的使用详解

gateway之断言的使用详解

文章目录 gateway产生的背景&#xff0c;为什么要是用gateway什么是网关gateway 带来的好处功能特征gateway在项目中使用的依赖 什么是断言断言分类内置自定义示例 断言和过滤器的不同 gateway产生的背景&#xff0c;为什么要是用gateway 一个系统会被拆分为多个微服务&#x…
阅读更多...
python处理CSV文件

python处理CSV文件

CSV库还有其他处理CSV的方法&#xff0c;这里只是介绍几个常用的&#xff0c;后面如果用到别的会进行更新 目录 1 生成一个新的csv文件&#xff0c;并向其中写一点东西 2 单纯往里面写几行 3 读取csv文件 1 生成一个新的csv文件&#xff0c;并向其中写一点东西 import…
阅读更多...
Java：关于mybatis框架mapper.xml编写小于号＜的问题

Java：关于mybatis框架mapper.xml编写小于号＜的问题

目录 方案一&#xff1a;转义字符方案二&#xff1a;原样字符总结参考文章 xml中小于号< 和 小于等于< 不能直接使用 select * from tb_user where age < #{user.age};方案一&#xff1a;转义字符 使用转义字符 含义符号转义字符小于<<大于>> 示例 s…
阅读更多...
最新文章

Copyright @ 2022~2023