跨境电商自建站平台,快速做网站公司报价,企业邮箱登录入口免费,公众号搭建2023 年是我最后一次参加 高教杯大学生数学建模竞赛 以后不会再参加了#xff08;大四参加意义不太#xff0c;研究生有研究生的数学建模大赛#xff09; 很遗憾 由于各种原因 我们没有能够完成赛题2022 年 美赛 2022年 Mathor Cup 2022 年国赛 2022 亚太杯 2023年 美赛 202…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, Lambda1, qv lambda t, x, y:0,X_start0, X_end1,Y_start0, Y_end1,T_start0, T_end1,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][Xi1][Yi])\2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]\ry * (self.results[Ti-1][Xi][Yi-1]self.results[Ti-1][Xi][Yi1])-\self.results[Ti-2][Xi][Yi]return self.resultsdef show_wave(self):fig plt.figure(figsize(8,8))ax1 fig.add_subplot(111, projection3d)x, y np.meshgrid(self.X, self.Y)for time in range(len(self.T)):ax1.plot_surface(x, y, self.results[time, :, :],rstride2, cstride2, cmaprainbow)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_start0, X_end1,Y_start0, Y_end1,T_start0, T_end1.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 test1: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][Xi1][Yi])\2* (1-rx-ry) *self.results[Ti-1][Xi][Yi]\ry * (self.results[Ti-1][Xi][Yi-1]self.results[Ti-1][Xi][Yi1])-\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, projection3d)ax2 fig.add_subplot(122, projection3d)x, y np.meshgrid(self.X, self.Y)for time in range(len(self.T)):ax1.plot_surface(x, y, self.results[time, :, :],rstride2, cstride2, cmaprainbow)ax2.plot_wireframe(x, y, self.results[time, :, :],rstride2, cstride2, linewidth1, cmaprainbow)plt.title(time: str(self.T[time]))plt.pause(0.01)plt.cla()def test_fun(x, y, t):return 1c waveequation2D(ctest_fun)
u c.Upwind3P()
c.show_wave()所以一般建议使用Matlab 中的pdetool 求解这是我写的两篇博客可以参考
典型的偏微分方程数值解法
二维Poisson方程五点差分格式与Python实现 一定要仔细阅读有关书籍否则很难做出比较好的成果
参考书目
偏微分方程的数值解法第三版 陆金甫特殊函数概论 王竹溪数学物理方法与仿真第三版 杨华军科学计算中的偏微分方程有限差分法 张文生 优化算法
这是是数学建模的核心优化求解一般而言很多问题可以做凸优化 但是数学建模中的优化问题往往不能...所以要么二分法全局搜索要么智能求解吧最近会写一些相关的博客就不具体演示了。也可以去我的博客下面找...... 写完博客心情舒畅不少
