高教杯数学建模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数值积分方法

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

中心差分式

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

数值偏微分

• 这也是常考的问题
• 但是数值偏微分方程的求解极其得复杂，稳定性特别难以分析，我举一个 用迎风法求解抛物型方程的例子

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实现

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

参考书目

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

优化算法

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

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