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

    • 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数值积分方法

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

    中心差分式

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

    数值偏微分

    • 这也是常考的问题
    • 但是数值偏微分方程的求解极其得复杂，稳定性特别难以分析，我举一个 用迎风法求解抛物型方程的例子
    1. import numpy as np
    2. import matplotlib.pyplot as plt
    3.  
    4. #\frac{\partial u}{\partial t} = \lambda(
    5. #\frac{\partial^2 u}{\partial x^2} +
    6. #\frac{\partial^2 u}{\partial y^2}
    7. #) + q_v
    8.  
    9.  
    10. class projequation2D():
    11.     def __init__(self, Lambda=1, qv = lambda t, x, y:0,
    12.                               X_start=0, X_end=1,
    13.                               Y_start=0, Y_end=1,
    14.                               T_start=0, T_end=1,
    15.                               dx = 0.05,
    16.                               dy = 0.05,
    17.                               dt = 0.01):
    18.         
    19.         #c = lambda x,y,t:(?)
    20.         self.Lambda = Lambda
    21.         self.qv = qv
    22.  
    23.         self.dx = dx
    24.         self.dy = dy
    25.         self.dt = dt
    26.  
    27.         self.X_start = X_start
    28.         self.Y_start = Y_start
    29.         self.T_start = T_start
    30.  
    31.         self.X_end = X_end
    32.         self.Y_end = Y_end
    33.         self.T_end = T_end
    34.  
    35.         self.X = np.arange(X_start, X_end, dx)
    36.         self.Y = np.arange(Y_start, Y_end, dy)
    37.         self.T = np.arange(T_start, T_end, dt)
    38.  
    39.         self.results = np.zeros((len(self.T), len(self.X), len(self.Y)))
    40.         self.startresults()
    41.  
    42.  
    43.     def initcondition_0(self):
    44.         X_num = len(self.X)
    45.         Y_num = len(self.Y)
    46.         test_fun = lambda x,y: np.sin(5*x)*np.cos(5*y)
    47.         for ix in range(X_num):
    48.             for iy in range(Y_num):
    49.                 self.results[0][ix][iy] == test_fun(self.X[ix], self.Y[iy])
    50.  
    51.     def boundarycondition_left(self):
    52.         Y_num = len(self.Y)
    53.         T_num = len(self.T)
    54.  
    55.         for it in range(T_num):
    56.             for iy in range(Y_num):
    57.                 self.results[it][0][iy] = 25
    58.  
    59.  
    60.         
    61.     def boundarycondition_right(self): 
    62.         Y_num = len(self.Y)
    63.         T_num = len(self.T)
    64.  
    65.         for it in range(T_num):
    66.             for iy in range(Y_num):
    67.                 self.results[it][-1][iy] = 25
    68.                 
    69.     def boundarycondition_up(self):
    70.         X_num = len(self.X)
    71.         T_num = len(self.T)
    72.  
    73.         for it in range(T_num):
    74.             for iy in range(X_num):
    75.                 self.results[it][iy][-1] = 25
    76.                 
    77.     def boundarycondition_down(self):
    78.         X_num = len(self.X)
    79.         T_num = len(self.T)
    80.  
    81.         for it in range(T_num):
    82.             for iy in range(X_num):
    83.                 self.results[it][iy][0] = 25
    84.         
    85.     def startresults(self):
    86.         self.initcondition_0()
    87.         self.boundarycondition_left()
    88.         self.boundarycondition_right()
    89.         self.boundarycondition_up()
    90.         self.boundarycondition_down()
    91.  
    92.     def Upwind3P(self):
    93.         for Ti in range(1, len(self.T)):
    94.             for Xi in range(1, len(self.X)-1):
    95.                 for Yi in range(1, len(self.Y)-1):
    96.                     rx = self.Lambda * self.dt/(self.dx)**2
    97.                     ry = self.Lambda * self.dt/(self.dy)**2
    98.                     self.results[Ti][Xi][Yi] = \
    99.                         rx * (self.results[Ti-1][Xi-1][Yi]+self.results[Ti-1][Xi+1][Yi])+\
    100.                         2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]+\
    101.                         ry * (self.results[Ti-1][Xi][Yi-1]+self.results[Ti-1][Xi][Yi+1])-\
    102.                         self.results[Ti-2][Xi][Yi]
    103.         return self.results
    104.     
    105.     def show_wave(self):
    106.         fig = plt.figure(figsize=(8,8))
    107.         ax1 = fig.add_subplot(111, projection='3d')
    108.         x, y = np.meshgrid(self.X, self.Y)
    109.         for time in range(len(self.T)):
    110.             ax1.plot_surface(x, y, self.results[time, :, :],
    111.                              rstride=2, cstride=2, cmap='rainbow')
    112.             plt.title("time:" + str(self.T[time]))
    113.             plt.pause(0.5)
    114.             plt.cla()
    115.  
    116. c = projequation2D()
    117. u = c.Upwind3P()
    118. c.show_wave()
    • 再来一个解波动方程的例子

    1. import numpy as np
    2. import matplotlib.pyplot as plt
    3.  
    4. #\frac{\partial^2 u}{\partial t^2} = c^2 (
    5. #\frac{\partial^2 u}{\partial x^2} +
    6. #\frac{\partial^2 u}{\partial y^2}
    7. #)
    8.  
    9.  
    10. class waveequation2D():
    11.     def __init__(self, c,
    12.                               X_start=0, X_end=1,
    13.                               Y_start=0, Y_end=1,
    14.                               T_start=0, T_end=1.01,
    15.                               dx = 0.01,
    16.                               dy = 0.01,
    17.                               dt = 0.01):
    18.         
    19.         #c = lambda x,y,t:(?)
    20.         self.c = c
    21.  
    22.         self.dx = dx
    23.         self.dy = dy
    24.         self.dt = dt
    25.  
    26.         self.X_start = X_start
    27.         self.Y_start = Y_start
    28.         self.T_start = T_start
    29.  
    30.         self.X_end = X_end
    31.         self.Y_end = Y_end
    32.         self.T_end = T_end
    33.  
    34.         self.X = np.arange(X_start, X_end, dx)
    35.         self.Y = np.arange(Y_start, Y_end, dy)
    36.         self.T = np.arange(T_start, T_end, dt)
    37.  
    38.         self.results = np.zeros((len(self.T), len(self.X), len(self.Y)))
    39.         self.startresults()
    40.  
    41.  
    42.     def initcondition_0(self, x, y):
    43.         return np.sin(10*x)
    44.  
    45.     def initcondition_1(self, x, y):
    46.         return np.sin(10*x)
    47.     
    48.     def boundarycondition_left(self, t, x, y):
    49.         return (np.sin(10*y))[0]
    50.  
    51.     def boundarycondition_right(self, t, x, y): 
    52.         return (np.sin(10*y))[0]
    53.  
    54.     def boundarycondition_up(self, t, x, y):
    55.         return (np.sin(10*x))[0]
    56.  
    57.     def boundarycondition_down(self, t, x, y):
    58.         return (np.sin(10*x))[0]
    59.  
    60.         
    61.     def startresults(self):
    62.         x, y = np.meshgrid(self.X, self.Y)
    63.         self.results[0] = self.initcondition_0(x, y)
    64.         self.results[1] = self.initcondition_1(x, y)
    65.  
    66.         t, x, y = np.meshgrid(self.T, self.X, self.Y)
    67.         self.results[:, 0, :]  = self.boundarycondition_left(t, self.X_start, y)
    68.         self.results[:, -1, :] = self.boundarycondition_right(t, self.X_end, y)
    69.         self.results[:, :, -1] = self.boundarycondition_up(t, x, self.Y_end)
    70.         self.results[:, :, 0]  = self.boundarycondition_down(t, x, self.Y_start)
    71.         
    72.         
    73.     def test_stability(self, x, y, t):
    74.         test = 4*self.dt**2*self.c(x, y, t)**2/(self.dx**2 + self.dy**2)
    75.         if test<=1:
    76.             return True
    77.         else:
    78.             print("unstable in x:", x, "y:", y, "t:", t)
    79.             return False
    80.  
    81.     def Upwind3P(self):
    82.         for Ti in range(2, len(self.T)):
    83.             for Xi in range(1, len(self.X)-1):
    84.                 for Yi in range(1, len(self.Y)-1):
    85.                     rx = self.c(self.X[Xi], self.Y[Yi], self.T[Ti])**2 * self.dt**2/self.dx**2
    86.                     ry = self.c(self.X[Xi], self.Y[Yi], self.T[Ti])**2 * self.dt**2/self.dy**2
    87.                     self.results[Ti][Xi][Yi] = \
    88.                         rx * (self.results[Ti-1][Xi-1][Yi]+self.results[Ti-1][Xi+1][Yi])+\
    89.                         2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]+\
    90.                         ry * (self.results[Ti-1][Xi][Yi-1]+self.results[Ti-1][Xi][Yi+1])-\
    91.                         self.results[Ti-2][Xi][Yi]
    92. ##                    self.test_stability(self.X[Xi],self.Y[Yi],self.T[Ti])                    
    93.         return self.results
    94.  
    95.     def show_wave(self):
    96.         fig = plt.figure(figsize=(8,12))
    97.         ax1 = fig.add_subplot(121, projection='3d')
    98.         ax2 = fig.add_subplot(122, projection='3d')
    99.  
    100.  
    101.         x, y = np.meshgrid(self.X, self.Y)
    102.         for time in range(len(self.T)):
    103.             ax1.plot_surface(x, y, self.results[time, :, :],
    104.                              rstride=2, cstride=2, cmap='rainbow')
    105.             ax2.plot_wireframe(x, y, self.results[time, :, :],
    106.                              rstride=2, cstride=2, linewidth=1, cmap='rainbow')
    107.             plt.title("time:" + str(self.T[time]))
    108.             plt.pause(0.01)
    109.             plt.cla()
    110.           
    111. def test_fun(x, y, t):
    112.     return 1
    113.  
    114.  
    115. c = waveequation2D(c=test_fun)
    116. u = c.Upwind3P()
    117. c.show_wave()

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

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

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

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

    参考书目

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

    优化算法

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

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

  • 相关阅读:
    利用Python分析金融交易中的滚动Z值
    五、kuternetes Pod介绍与配置
    【计算机组成原理】习题（一）—— 指令系统
    2022烟台海洋牧场优势产品宣传推介会在杭州成功举办
    基于JAVA校园二手书交易系统计算机毕业设计源码+数据库+lw文档+系统+部署
    判断链表是否是环形链表
    MySQL性能调优
    leetcode 刷题 log day 46（背包总结篇
    python爬虫，多线程与生产者消费者模式
    LeetCode题解—260.只出现一次的数字Ⅲ
  • 原文地址：https://blog.csdn.net/Chandler_river/article/details/133053545