高教杯数学建模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()