🎓 本文收录于 《滚雪球学数学建模》专栏，本专栏专注于从零基础出发，带你循序渐进掌握数学建模的核心方法与思维路径。通过“滚雪球式”学习法，从最简单的生活化案例到经典建模工具，再到省赛、国赛真题实战演练，逐步积累知识与技能，帮助你搭建完整的建模体系。无论你是刚接触建模的新生，还是希望提升科研与竞赛能力的学者与工程师，本专栏都将为你提供一条清晰、高效、可实操的建模成长路线，助力你快速进阶，突破难点，迈向数学建模竞赛与实际应用的更高峰 🚀！
全文目录：
🔧 前言
哎呀呀，亲爱的建模探险家们，咱们的工具箱挖掘之旅又要掀起新一轮高潮了！在上个部分，我像个没完没了的探险队长，领着你们在概率统计的迷宫里转悠，那些分布的曲线和期望的计算，是不是让你们觉得随机世界突然间变得那么亲切而有趣？哈哈，我自己写着写着就回想起那些疯狂的日子：有次凌晨三点，盯着屏幕上的方差图，突然灵光一闪，调整了个参数，模型从一塌糊涂变成完美无缺，那种“啊哈”时刻的喜悦，简直想跳起来拥抱电脑！当然，也少不了那些坑人的瞬间，比如估计偏差太大，导致整个预测翻车，我沮丧得差点砸键盘。
但也正是这些起起落落，让我从一个数学小白成长为动态建模的“老司机”。现在，继续咱们的动态模型与微分方程这一章，我得认真响应你们的反馈，内容再爆炸式丰富，字数堆成小山，让你们读得过瘾到停不下来。深度上，我不光停留在公式表面，还挖到数值方法的误差分析、混沌系统的敏感依赖、甚至量子力学中的薛定谔方程类比；
广度上，从细胞内分子扩散的小尺度，到全球气候系统的宏大循环，多维度拉扯开来，跨界到神经科学、股票市场波动、甚至用方程模拟人际关系的“情感动态”（没错，我曾经试过用一个简单的阻尼振荡模型分析友情起伏，结果意外地准，哈哈，但这只是娱乐，别太认真哦！）。
情感上说，这些动态工具不只是冰冷的数学符号，它们是我人生中的“生命之河”，帮我捕捉变化的本质，解决过无数现实难题，从帮朋友优化生意流程到自己反思人生的起伏曲线。哎，别磨蹭了，我会多洒私货：更多我的囧事、心得小窍门、灵感闪现点，保证读起来像听老友拉家常，笑声不断，启发无限。来吧，深吸一口气，准备让时间在你的模型中“活蹦乱跳”——你敢打赌，这些方程不会让你重新审视自己的人生“轨迹曲线”吗？哈哈，一起冲刺吧！
正文
🚀 动态模型与微分方程
哇哦，这一大块绝对是我工具箱里的“时间引擎”！动态模型与微分方程，不是那种让人脑仁疼的抽象理论，而是帮你模拟现实世界如何随时间演变的超级魔法杖。生活中处处是动态变化：心跳的节律、河流的奔腾、股市的涨跌、疫情的扩散——没有这些工具，你的建模就卡在“定格瞬间”，缺少了那份活生生的脉动和故事性。我回想自己第一次真正爱上微分方程，是在大学的一个项目中，模拟一个弹簧的振动，那解出来的正弦曲线让我激动得夜不能寐，感觉自己像发现了宇宙的秘密节奏！但也出过天大的糗事：有次建模一个化学反应系统，忽略了初始条件的敏感性，方程解发散成天文数字，我调试了整整两天，累到眼睛红肿，但最终搞定后，那种征服感的满足，啧啧，值了！
这次，我把内容升级到极致：基础部分多加历史演变、多解法对比、多技巧扩展；案例部分不只原有的两个，我多加了三个（捕食者-猎物、弹簧振动、化学反应、额外加传染病传播和股票波动），每个案例都配独立MATLAB代码，一行行详细解析，确保新手也能上手。深度上，层层剖析解析解的局限性、数值解的误差控制、平衡点的李雅普诺夫稳定性、参数敏感性分析，甚至触及分岔理论和混沌吸引子；广度上，横跨三十多个领域，从微观生物细胞分裂到宏观宇宙黑洞蒸发模拟。哈哈，保证读完后你会感慨万千：变化才是宇宙的永恒主题，掌握动态模型，你就是时间的“主宰者”！准备好让你的建模世界“动起来”了吗？哎，人生不也像一条复杂的微分曲线，充满未知的转折和惊喜？
1. 常微分方程与变化建模
嘿嘿，先来彻底拆解常微分方程（ODE）这个“动态核心”吧！常微分方程本质上是描述一个或多个变量对独立变量（通常是时间t）变化率的数学方程，通用形式包括一阶： d y d t = f ( t , y ) \frac{dy}{dt} = f(t, y) dtdy=f(t,y)，或更高阶如二阶： d 2 y d t 2 = g ( t , y , d y d t ) \frac{d^2 y}{dt^2} = g(t, y, \frac{dy}{dt}) dt2d2y=g(t,y,dtdy)。它在建模变化时超级强大，因为它把“率”作为主角，让你从静态描述跃升到动态预测。简单来说，像一个物体的自由落体运动 d v d t = g \frac{dv}{dt} = g dtdv=g（v速度，g重力加速度），积分得 v = g t + v 0 v = gt + v_0 v=gt+v0，再积分位置 s = 1 2 g t 2 + v 0 t + s 0 s = \frac{1}{2} gt^2 + v_0 t + s_0 s=21gt2+v0t+s0，预测抛物线轨迹。为什么我对ODE情有独钟？因为它让我从死记硬背公式转向真正理解“过程的流动”，像一个艺术家用画笔捕捉光影的变幻。我记得入门阶段，解一个简单的指数衰减方程 d y d t = − k y \frac{dy}{dt} = -k y dtdy=−ky，得到 y = y 0 e − k t y = y_0 e^{-kt} y=y0e−kt，那优雅的解析形式让我脑洞大开，兴奋得立刻应用到模拟咖啡冷却上！但也闹过大笑话：有次建模一个非线性系统 d y d t = y ( 1 − y ) + sin t \frac{dy}{dt} = y (1 - y) + \sin t dtdy=y(1−y)+sint，盲目用线性近似，结果预测偏差巨大，我在报告会上被导师点名批评，尴尬得想找地缝钻，但那次失败让我深入学习了摄动方法，从此模型更robust。
深度上，咱们多挖几层泥土，先分类聊聊：按阶数，一阶最基础，如分离变量法 ∫ d y f ( y ) = ∫ g ( t ) d t \int \frac{dy}{f(y)} = \int g(t) dt ∫f(y)dy=∫g(t)dt解出显式形式；二阶常见于物理，如简谐振子 d 2 x d t 2 + ω 2 x = 0 \frac{d^2 x}{dt^2} + \omega^2 x = 0 dt2d2x+ω2x=0，特征方程 r 2 + ω 2 = 0 r^2 + \omega^2 = 0 r2+ω2=0得 r = ± i ω r = \pm i \omega r=±iω，解 x = A cos ω t + B sin ω t x = A \cos \omega t + B \sin \omega t x=Acosωt+Bsinωt；更高阶如三阶系统在流体力学中 d 3 y d t 3 + a d 2 y d t 2 + b d y d t + c y = 0 \frac{d^3 y}{dt^3} + a \frac{d^2 y}{dt^2} + b \frac{dy}{dt} + c y = 0 dt3d3y+adt2d2y+bdtdy+cy=0，用伴随多项式求根。线性与非线性大不同：线性ODE叠加原理适用 y = y h + y p y = y_h + y_p y=yh+yp（齐次 + 特解），如未定系数法求特解 y p = A e r t y_p = A e^{rt} yp=Aert；非线性如Van der Pol振荡器 d 2 x d t 2 − μ ( 1 − x 2 ) d x d t + x = 0 \frac{d^2 x}{dt^2} - \mu (1 - x^2) \frac{dx}{dt} + x = 0 dt2d2x−μ(1−x2)dtdx+x=0，引入限周期，深度上用相平面分析 x ˙ = v , v ˙ = μ ( 1 − x 2 ) v − x \dot x = v, \dot v = \mu (1 - x^2) v - x x˙=v,v˙=μ(1−x2)v−x，挖到吸引子概念，甚至链接到混沌理论的奇异吸引子。我特别爱非线性，因为它捕捉现实的复杂性，比如Duffing方程 d 2 x d t 2 + δ d x d t + α x + β x 3 = γ cos ω t \frac{d^2 x}{dt^2} + \delta \frac{dx}{dt} + \alpha x + \beta x^3 = \gamma \cos \omega t dt2d2x+δdtdx+αx+βx3=γcosωt，分岔行为丰富，能模拟从有序到混沌的转变，我用Poincaré截面分析过，深度让我感慨自然的不可预测性。
解法扩展多点：解析解有限，如Bernoulli方程 d y d t + P y = Q y n \frac{dy}{dt} + P y = Q y^n dtdy+Py=Qyn，置换 v = y 1 − n v = y^{1-n} v=y1−n转线性；更多靠数值，如Euler方法基础但截断误差O(h)，改进Heun k 1 = f ( t , y ) , k 2 = f ( t + h , y + h k 1 ) , y n + 1 = y n + h 2 ( k 1 + k 2 ) O ( h 2 ) k1 = f(t,y), k2 = f(t+h, y+h k1), y_{n+1} = y_n + \frac{h}{2} (k1 + k2) O(h^2) k1=f(t,y),k2=f(t+h,y+hk1),yn+1=yn+2h(k1+k2)O(h2)；金字塔是Runge-Kutta家族，RK4高精度O(h^4)，公式 k 1 = h f ( t , y ) , k 2 = h f ( t + h / 2 , y + k 1 / 2 ) , k 3 = h f ( t + h / 2 , y + k 2 / 2 ) , k 4 = h f ( t + h , y + k 3 ) , y n + 1 = y n + ( k 1 + 2 k 2 + 2 k 3 + k 4 ) / 6 k1 = h f(t,y), k2 = h f(t+h/2, y+k1/2), k3 = h f(t+h/2, y+k2/2), k4 = h f(t+h, y+k3), y_{n+1} = y_n + (k1 + 2k2 + 2k3 + k4)/6 k1=hf(t,y),k2=hf(t+h/2,y+k1/2),k3=hf(t+h/2,y+k2/2),k4=hf(t+h,y+k3),yn+1=yn+(k1+2k2+2k3+k4)/6;
我常用在刚性问题；刚性ODE如化学快速反应，用隐式方法如backward Euler y n + 1 = y n + h f ( t n + 1 , y n + 1 ) y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}) yn+1=yn+hf(tn+1,yn+1)，牛顿迭代求解稳定性好。稳定性分析深度挖：平衡点 f ( y ∗ ) = 0 f(y^*) = 0 f(y∗)=0，线性化雅可比矩阵 J = ∂ f ∂ y ∣ y ∗ J = \frac{\partial f}{\partial y} |_{y^*} J=∂y∂f∣y∗，特征值Re(\lambda)<0收敛；李雅普诺夫函数 V ( y ) > 0 , V ˙ < 0 V(y) >0, \dot V <0 V(y)>0,V˙<0证全局稳定，我用 V = ( y − y ∗ ) 2 / 2 V = (y - y^*)^2 /2 V=(y−y∗)2/2简单判。敏感性分析 S = ∂ y ∂ p / ( y / p ) S = \frac{\partial y}{\partial p} / (y / p) S=∂p∂y/(y/p)评估参数影响，用有限差分或伴随方法计算，深度上帮优化设计。
历史典故多添几个生动点：牛顿莱布尼茨发明微积分就是为动态问题，如牛顿运动定律 m d 2 r d t 2 = F m \frac{d^2 \mathbf{r}}{dt^2} = \mathbf{F} mdt2d2r=F奠基力学；Euler发展多体问题，Lagrange引入变分法 δ ∫ L d t = 0 \delta \int L dt = 0 δ∫Ldt=0（L动能-势能），现代控制理论基础；Poincaré开创拓扑动力学，Lorenz 1963年发现混沌 x ˙ = σ ( y − x ) , y ˙ = x ( ρ − z ) − y , z ˙ = x y − β z \dot x = \sigma (y - x), \dot y = x (\rho - z) - y, \dot z = x y - \beta z x˙=σ(y−x),y˙=x(ρ−z)−y,z˙=xy−βz，蝴蝶效应让我模拟天气时，深度敬畏初始条件的敏感依赖。广度上拉得老长了：在生物学，FitzHugh-Nagumo神经元 v ˙ = v − v 3 / 3 − w + I , w ˙ = ϵ ( v + a − b w ) \dot v = v - v^3/3 - w + I, \dot w = \epsilon (v + a - b w) v˙=v−v3/3−w+I,w˙=ϵ(v+a−bw)模拟脉冲；物理，阻尼驱动 x ¨ + 2 ζ ω x ˙ + ω 2 x = F cos ω d t \ddot x + 2 \zeta \omega \dot x + \omega^2 x = F \cos \omega_d t x¨+2ζωx˙+ω2x=Fcosωdt，共振放大 A = F / ( ω 2 − ω d 2 ) 2 + ( 2 ζ ω ω d ) 2 A = F / \sqrt{( \omega^2 - \omega_d^2 )^2 + (2 \zeta \omega \omega_d)^2} A=F/(ω2−ωd2)2+(2ζωωd)2；经济学，哈罗德-多马增长 k ˙ = s f ( k ) − ( n + δ ) k \dot k = s f(k) - (n + \delta) k k˙=sf(k)−(n+δ)k，稳态 k ∗ = ( s / ( n + δ ) ) 1 / ( 1 − α ) k^* = (s / (n + \delta))^{1/(1-\alpha)} k∗=(s/(n+δ))1/(1−α)，我用模拟发展中国家经济，帮政策分析，情感那种用数学推动社会进步的感觉，太暖心了！工程，PID控制 u = K p e + K i ∫ e d t + K d e ˙ u = K_p e + K_i \int e dt + K_d \dot e u=Kpe+Ki∫edt+Kde˙基于ODE稳定性；医学，药物代谢 C ˙ = − k C + D δ ( t − t d ) \dot C = -k C + D \delta (t - t_d) C˙=−kC+Dδ(t−td)多剂量，预测血浓峰谷。
变化建模技巧多积累：先识别关键率，如速度
x
˙
=
v
\dot x = v
x˙=v、加速度
v
˙
=
F
/
m
\dot v = F/m
v˙=F/m；假设简化如小角度近似
sin
θ
≈
θ
\sin \theta \approx \theta
sinθ≈θ；数据验证用卡尔曼滤波融合噪声
d
o
t
x
^
=
A
x
^
+
K
(
z
−
C
x
^
)
dot \hat x = A \hat x + K (z - C \hat x)
dotx^=Ax^+K(z−Cx^)；参数估计用遗传算法优
m
i
n
J
=
∫
(
y
m
o
d
e
l
−
y
d
a
t
a
)
2
d
t
min J = \int (y_{model} - y_{data})^2 dt
minJ=∫(ymodel−ydata)2dt。实用工具扩展：MATLAB ode45自适应RK，
o
d
e
s
e
t
(
′
A
b
s
T
o
l
′
,
1
e
−
8
,
′
R
e
l
T
o
l
′
,
1
e
−
6
)
odeset('AbsTol',1e-8,'RelTol',1e-6)
odeset(′AbsTol′,1e−8,′RelTol′,1e−6) 控误差；Python scipy.integrate.solve_ivp多方法选。误区扒皮多层：初始条件小扰大影响，混沌系统需Lorenz地图；步长固定忽略局部误差，用嵌入RK如ode45自动调；稳定性盲判忽略鞍点，用分岔图
μ
\mu
μ变看 Hopf 分岔
λ
=
α
±
i
β
\lambda = \alpha \pm i \beta
λ=α±iβ生振荡。个人情感大吐槽：ODE让我悟到人生无常——像非线性系统，小选择大结局，有次模拟职业路径
c
˙
=
r
c
(
1
−
c
/
K
)
−
d
\dot c = r c (1 - c/K) - d
c˙=rc(1−c/K)−d，帮自己规划，转行成功，那自省的温暖，终身难忘。
更多故事：用振动方程 θ ¨ + ( g / l ) θ = 0 \ddot \theta + (g/l) \theta = 0 θ¨+(g/l)θ=0建模钟摆帮侄子物理作业，他考满分，叔侄欢呼！多维度拓展：AI中连续时间RNN用神经ODE h ˙ = f ( h , t , θ ) \dot h = f(h, t, \theta) h˙=f(h,t,θ)学序列；伦理角度，动态模拟基因编辑 g ˙ = m g ( 1 − g / K ) − d \dot g = m g (1 - g/K) - d g˙=mg(1−g/K)−d评估风险防失控；厨艺实践，用加热方程 T ˙ = h ( T h − T ) + q \dot T = h (T_h - T) + q T˙=h(Th−T)+q优烤饼时间，避免焦黑；社会学，用SIR扩展 S ˙ = − β S I + ν R \dot S = -\beta S I + \nu R S˙=−βSI+νR模拟意见传播，帮营销策略。广度再拓宽：神经科学，Hodgkin-Huxley V ˙ = ( I − g N a m 3 h ( V − E N a ) − g K n 4 ( V − E K ) − g l ( V − E l ) ) / C \dot V = (I - g_{Na} m^3 h (V - E_{Na}) - g_K n^4 (V - E_K) - g_l (V - E_l))/C V˙=(I−gNam3h(V−ENa)−gKn4(V−EK)−gl(V−El))/C离子通道；地质，热扩散 T ˙ = κ ∇ 2 T \dot T = \kappa \nabla^2 T T˙=κ∇2T地震预测；艺术，用Lissajous曲线 x = A sin ( a t + ϕ ) , y = B sin b t x = A \sin (a t + \phi), y = B \sin b t x=Asin(at+ϕ),y=Bsinbt生成动态图案，我试过做屏保，创意爆棚！哎呀，ODE的世界无穷尽，多内容让你眼界大开，但记住秘诀：从简单一阶起步，多数值模拟多稳定性检查，你会渐渐爱上这个“华丽转身”的艺术！
2. 案例：人口增长、咖啡降温、捕食者-猎物、弹簧振动、化学反应、传染病传播、股票波动
哈哈，理论吃饱了，来多案例的实战狂欢吧！我不光原有两个，多加五个，总七个，每个案例深度剖析背景、方程推导、参数设置、结果解读、扩展讨论、我的亲身故事，再配独立MATLAB代码，一行行解析，确保你能复制跑起来。内容丰富如盛宴，从生态到金融，从物理到生物，多维度展示动态魅力。
代码1：人口增长logistic MATLAB版
第一个：人口增长。经典Malthus模型 P ˙ = r P \dot P = r P P˙=rP指数增长，但现实有限制，用logistic P ˙ = r P ( 1 − P / K ) \dot P = r P (1 - P/K) P˙=rP(1−P/K)，S形曲线收敛到承载力K。我模拟一个小镇人口，r=0.018, K=2e5, P0=5e4，预测10年内峰值，帮当地政府规划基础设施。但忽略移民流动，曲线偏高，扩展加迁移项 + m − e P + m - e P +m−eP准了。情感上，建模人口让我感慨生命的韧性——增长有限，需可持续发展。
% 人口增长logistic模型 - 小明老师实战版
% 功能: 数值解ODE，画时间曲线，分析稳态
clear; clc; close all;
r = 0.018; K = 2e5; P0 = 5e4; tspan = [0 50]; % 年单位
pop_ode = @(t,P) r * P * (1 - P/K); % 定义变化率 dP/dt
[t, P] = ode45(pop_ode, tspan, P0); % 用RK方法数值求解
figure; plot(t, P/1e3, 'b-', 'LineWidth', 1.5); % 画曲线，人口单位千
xlabel('时间 (年)'); ylabel('人口 (千)'); title('人口增长动态 - Logistic模型'); grid on;
hold on; plot(tspan, K/1e3 * ones(size(tspan)), 'r--'); legend('P(t)', '承载力K'); % 加水平线示K
% 扩展: 加扰动模拟现实不确定
P_perturb = P + 1e3 * randn(size(P)); plot(t, P_perturb/1e3, 'g:'); legend('P(t)', 'K', '加噪模拟'); % 加随机噪
解析： 清变量初始化参数现实值；ode函数定义logistic率；ode45自动步长解轨迹；画图加承载线直观收敛；扩展加randn噪模拟不确定，如出生率波动。跑这个，调r见爆炸 vs 稳定，超级实用预测城市压力！
代码2：咖啡降温Newton加拟合与扩展
第二个：咖啡降温。Newton冷却定律 T ˙ = − k ( T − T a ) \dot T = -k (T - T_a) T˙=−k(T−Ta)，T温度，Ta环境，k冷却常数。解析解 T = T a + ( T 0 − T a ) e − k t T = T_a + (T_0 - T_a) e^{-kt} T=Ta+(T0−Ta)e−kt。我实际实验杯子咖啡，测T0=95°C, Ta=20°C, k=0.05，预测最佳饮温60°C时间约15分。但忽略对流风速，k变小，扩展加Stefan-Boltzmann辐射 − ϵ σ ( T 4 − T a 4 ) - \epsilon \sigma (T^4 - T_a^4) −ϵσ(T4−Ta4)准辐射效应。广度类似热交换工程，情感让我珍惜“热”时刻——咖啡凉得快，人生机会稍纵即逝。
% 咖啡降温Newton模型 - 加参数拟合 & 辐射扩展
clear; clc; close all;
k = 0.05; Ta = 20; T0 = 95; tspan = [0 30]; % 时间分钟
coffee_ode = @(t,T) -k * (T - Ta); % 基本冷却率
[t, T] = ode45(coffee_ode, tspan, T0); % 数值解
figure; plot(t, T, 'r-', 'LineWidth', 1.5); hold on; % 画基本曲线
plot(tspan, Ta * ones(size(tspan)), 'g--'); xlabel('时间 (分)'); ylabel('温度 (°C)');
title('咖啡降温动态 - Newton模型'); grid on; legend('T(t)', '环境Ta');
% 参数拟合 - 假设实验数据
t_data = [0 5 10 15 20 25 30]; T_data = [95 80 68 58 50 44 39]; % 测得数据
fit_fun = @(param, t) Ta + (T0 - Ta) * exp(-param(1)*t); % 解析形式拟合k
param0 = 0.03; % 初猜
k_fit = lsqcurvefit(fit_fun, param0, t_data, T_data); % 最小二乘拟合
T_fit = fit_fun(k_fit, t); plot(t, T_fit, 'b--'); legend('T(t)', 'Ta', '拟合曲线');
fprintf('拟合冷却常数k: %.4f\n', k_fit); % 输出k
% 扩展: 加辐射项 \epsilon \sigma (T^4 - Ta^4), 假设\epsilon\sigma = 1e-10
rad_ode = @(t,T) -k * (T - Ta) - 1e-10 * (T^4 - Ta^4); % 扩展率
[t_rad, T_rad] = ode45(rad_ode, tspan, T0); plot(t_rad, T_rad, 'm:'); legend('T(t)', 'Ta', '拟合', '加辐射');
解析： 参数设实验值；ode定义基本率；ode45解画加环境线；lsqcurvefit用解析拟合数据得k，画拟合验证；扩展rad_ode加T^4辐射项再解对比，见高温时辐射加速凉。数据从 t h e r m o m e t e r thermometer thermometer，跑调 ϵ \epsilon ϵ见影响，工程热设计神器！
代码3：捕食者-猎物Lotka-Volterra加相平面
第三个：捕食者-猎物。 L o t k a − V o l t e r r a Lotka-Volterra Lotka−Volterra模型 x ˙ = α x − β x y \dot x = \alpha x - \beta x y x˙=αx−βxy, y ˙ = − γ y + δ x y \dot y = -\gamma y + \delta x y y˙=−γy+δxy，x猎物y捕食者。参数 α = 0.08 , β = 0.002 , γ = 0.04 , δ = 0.001 , 初 x 0 = 200 , y 0 = 50 \alpha=0.08, \beta=0.002, \gamma=0.04, \delta=0.001, 初x0=200, y0=50 α=0.08,β=0.002,γ=0.04,δ=0.001,初x0=200,y0=50，周期振荡模拟生态循环。我用帮野生保护区，预测狼兔平衡，但忽略季节，振幅大，扩展加时间依赖 β ( t ) = β _ 0 ( 1 + sin 2 π t / 12 ) \beta(t) = \beta\_0 (1 + \sin 2\pi t /12) β(t)=β_0(1+sin2πt/12)准月变。广度生态学经典，情感如关系动态——互依共存，失衡灭顶。
% 捕食者-猎物Lotka-Volterra模型 - 加相平面 & 平衡点分析
clear; clc; close all;
alpha = 0.08; beta = 0.002; gamma = 0.04; delta = 0.001; x0 = 200; y0 = 50; tspan = [0 100]; % 时间任意单位
lv_ode = @(t,z) [alpha*z(1) - beta*z(1)*z(2); -gamma*z(2) + delta*z(1)*z(2)]; % z=[x y]向量ODE
[t, z] = ode45(lv_ode, tspan, [x0 y0]); % 数值解
figure(1); plot(t, z(:,1), 'g-', 'LineWidth', 1.5); hold on; plot(t, z(:,2), 'm-', 'LineWidth', 1.5);
xlabel('时间'); ylabel('种群数量'); title('捕食者-猎物动态'); grid on; legend('猎物x', '捕食者y');
% 相平面分析 - 轨迹闭合循环
figure(2); plot(z(:,1), z(:,2), 'k-', 'LineWidth', 2); xlabel('猎物x'); ylabel('捕食者y');
title('相平面 - 周期轨道'); grid on;
% 平衡点计算: x* = gamma/delta, y* = alpha/beta
x_star = gamma / delta; y_star = alpha / beta; hold on; plot(x_star, y_star, 'ro', 'MarkerSize', 10);
text(x_star+10, y_star+10, '平衡点'); % 标平衡
% 扩展: 加季节扰动 beta(t)
beta_t = @(t) beta * (1 + 0.2 * sin(2*pi*t/12)); % 月周期扰
lv_perturb = @(t,z) [alpha*z(1) - beta_t(t)*z(1)*z(2); -gamma*z(2) + delta*z(1)*z(2)];
[t_p, z_p] = ode45(lv_perturb, tspan, [x0 y0]); figure(3); plot(t_p, z_p(:,1), 'c:', t_p, z_p(:,2), 'y:');
title('加季节扰动动态'); legend('x扰', 'y扰');
解析：参数设生态现实； l v _ o d e lv\_ode lv_ode向量定义两率；ode45解画时序曲线见振荡；相平面 p l o t ( z ( : , 1 ) , z ( : , 2 ) ) plot(z(:,1),z(:,2)) plot(z(:,1),z(:,2))见闭环轨道；计算平衡点标红点；扩展 b e t a _ t beta\_t beta_t时间依再解对比，见扰动乱周期。跑调delta见灭绝风险，生态模拟利器！
代码4：弹簧振动阻尼谐振加驱动力
第四个：弹簧振动。阻尼谐振模型 x ¨ + 2 ζ ω x ˙ + ω 2 x = 0 \ddot x + 2 \zeta \omega \dot x + \omega^2 x = 0 x¨+2ζωx˙+ω2x=0， ω = k / m , ζ = c / ( 2 m k ) \omega = \sqrt{k/m}, \zeta = c/(2\sqrt{m k}) ω=k/m,ζ=c/(2mk)。参数 ω = 8 , ζ = 0.3 , x 0 = 0.5 , v 0 = 0 \omega=8, \zeta=0.3, x0=0.5, v0=0 ω=8,ζ=0.3,x0=0.5,v0=0，欠阻尼振荡衰减。我建模桥梁振动，预测共振风险，但忽略非线性，振幅过大，扩展加Duffing立方项 + β x 3 + \beta x^3 +βx3准硬弹簧。广度机械工程核心，情感如人生振荡——有起有落，阻尼帮稳。
% 弹簧振动阻尼谐振模型 - 二阶转一阶，加驱动力扩展
clear; clc; close all;
omega = 8; zeta = 0.3; x0 = 0.5; v0 = 0; tspan = [0 10]; % 时间秒
spring_ode = @(t,y) [y(2); -omega^2 * y(1) - 2*zeta*omega * y(2)]; % y=[x v], 基本无驱
[t, y] = ode45(spring_ode, tspan, [x0 v0]); % 解
figure(1); plot(t, y(:,1), 'c-', 'LineWidth', 1.5); xlabel('时间 (秒)'); ylabel('位移x (m)');
title('弹簧振动动态 - 欠阻尼'); grid on;
% 扩展: 加驱动力 F cos \omega_d t, \omega_d近 \omega 共振
F = 0.1; omega_d = 7.5; % 驱动幅 & 频
driven_ode = @(t,y) [y(2); -omega^2 * y(1) - 2*zeta*omega * y(2) + F * cos(omega_d * t)];
[t_d, y_d] = ode45(driven_ode, tspan, [x0 v0]);
figure(2); plot(t_d, y_d(:,1), 'k--'); title('加驱动力 - 近共振放大'); grid on;
% 稳定性分析: 特征值 r = -zeta \omega \pm \sqrt{(zeta \omega)^2 - \omega^2}
r1 = -zeta*omega + sqrt((zeta*omega)^2 - omega^2); r2 = -zeta*omega - sqrt((zeta*omega)^2 - omega^2);
fprintf('特征值: r1=%.2f + %.2fi, r2=%.2f + %.2fi\n', real(r1), imag(r1), real(r2), imag(r2)); % 复数示振荡
解析：参数设欠阻； s p r i n g _ o d e spring\_ode spring_ode转一阶向量； o d e 45 ode45 ode45解画位移衰振；扩展 d r i v e n _ o d e driven\_ode driven_ode加 c o s cos cos驱再解见放大；计算特征值示稳定（实部负） & 振荡（虚部）。跑调 z e t a zeta zeta见临界阻尼，工程防振必备！
代码5：化学反应二级加逆扩展
第五个：化学反应。二级反应A+B→C， A ˙ = − k A B , B ˙ = − k A B , C ˙ = k A B \dot A = -k A B, \dot B = -k A B, \dot C = k A B A˙=−kAB,B˙=−kAB,C˙=kAB。参数 k = 0.005 , A 0 = 2 , B 0 = 1.5 , C 0 = 0 k=0.005, A0=2, B0=1.5, C0=0 k=0.005,A0=2,B0=1.5,C0=0，浓度变化预测产量。我实验室模拟，拟合k帮优化催化，但忽略逆反应，平衡错，扩展加逆 + k r C + k_r C +krC准可逆。广度化工基础，情感反应如机遇——催化加速，逆转需防。
% 化学反应A+B→C模型 - 加逆反应扩展
clear; clc; close all;
k = 0.005; A0 = 2; B0 = 1.5; C0 = 0; tspan = [0 200]; % 时间任意
chem_ode = @(t,z) [-k*z(1)*z(2); -k*z(1)*z(2); k*z(1)*z(2)]; % z=[A B C]基本率
[t, z] = ode45(chem_ode, tspan, [A0 B0 C0]); % 解
figure; plot(t, z(:,1), 'r-', t, z(:,2), 'b-', t, z(:,3), 'g-', 'LineWidth', 1.5);
xlabel('时间'); ylabel('浓度 (mol/L)'); title('化学反应动态'); grid on; legend('A', 'B', 'C');
% 扩展: 加逆反应 k_r C → A+B, k_r=0.001
k_r = 0.001; rev_ode = @(t,z) [-k*z(1)*z(2) + k_r*z(3); -k*z(1)*z(2) + k_r*z(3); k*z(1)*z(2) - k_r*z(3)];
[t_r, z_r] = ode45(rev_ode, tspan, [A0 B0 C0]);
hold on; plot(t_r, z_r(:,1), 'r:', t_r, z_r(:,2), 'b:', t_r, z_r(:,3), 'g:'); legend('A', 'B', 'C', 'A逆', 'B逆', 'C逆');
解析：参数设初始浓度； c h e m _ o d e chem\_ode chem_ode向量三率；ode45解画变化（A B降 C升）；扩展 r e v _ o d e rev\_ode rev_ode加逆项再解对比，见平衡不零。跑调 k _ r k\_r k_r见可逆影响，化工优化神！
代码6：传染病SIR加免疫衰减。
第六个：传染病传播。SIR模型扩展 S ˙ = − β S I / N + ν R , I ˙ = β S I / N − γ I , R ˙ = γ I − ν R \dot S = - \beta S I / N + \nu R, \dot I = \beta S I / N - \gamma I, \dot R = \gamma I - \nu R S˙=−βSI/N+νR,I˙=βSI/N−γI,R˙=γI−νR，加恢复免疫衰减\nu。参数 β = 0.4 , γ = 0.1 , ν = 0.01 , N = 1 e 5 , S 0 = 0.99 N , I 0 = 0.01 N , R 0 = 0 \beta=0.4, \gamma=0.1, \nu=0.01, N=1e5, S0=0.99N, I0=0.01N, R0=0 β=0.4,γ=0.1,ν=0.01,N=1e5,S0=0.99N,I0=0.01N,R0=0，预测峰值。我疫情时用，帮社区预警，但忽略超级传播，I峰高，扩展加随机项准。广度流行病学经典，情感疫情如波澜——早建模早防控。
% 传染病传播SIR模型 - 加免疫衰减扩展
clear; clc; close all;
beta = 0.4; gamma = 0.1; nu = 0.01; N = 1e5; S0 = 0.99*N; I0 = 0.01*N; R0 = 0; tspan = [0 100]; % 天
sir_ode = @(t,z) [-beta*z(1)*z(2)/N; beta*z(1)*z(2)/N - gamma*z(2); gamma*z(2)]; % 基本SIR z=[S I R]
[t, z] = ode45(sir_ode, tspan, [S0 I0 R0]); % 解
figure; plot(t, z(:,1)/N, 'b-', t, z(:,2)/N, 'r-', t, z(:,3)/N, 'g-', 'LineWidth', 1.5);
xlabel('时间 (天)'); ylabel('比例'); title('传染病传播动态'); grid on; legend('S', 'I', 'R');
% 扩展: 加免疫衰减 \nu R → S
ext_ode = @(t,z) [-beta*z(1)*z(2)/N + nu*z(3); beta*z(1)*z(2)/N - gamma*z(2); gamma*z(2) - nu*z(3)];
[t_e, z_e] = ode45(ext_ode, tspan, [S0 I0 R0]);
hold on; plot(t_e, z_e(:,1)/N, 'b:', t_e, z_e(:,2)/N, 'r:', t_e, z_e(:,3)/N, 'g:'); legend('S', 'I', 'R', 'S衰', 'I衰', 'R衰');
解析： 参数设 R 0 = β / γ = 4 R0= \beta / \gamma =4 R0=β/γ=4高传； s i r _ o d e sir\_ode sir_ode向量三率；ode45解画比例曲线见I峰；扩展 e x t _ o d e ext\_ode ext_ode加 ν \nu ν循环再解对比，见免疫衰反复波。跑调 ν \nu ν见长期行为，疫情策略宝！
代码7：股票波动几何布朗加蒙特卡洛。
第七个：股票波动。随机游走扩展 S ˙ = μ S + σ S ξ \dot S = \mu S + \sigma S \xi S˙=μS+σSξ，几何布朗运动， μ \mu μ漂移, σ \sigma σ波动, ξ \xi ξ白噪。但连续化，用Ito积分近似。我模拟股价， m u = 0.1 mu=0.1 mu=0.1, σ = 0.2 , S 0 = 100 \sigma=0.2, S0=100 σ=0.2,S0=100，预测风险，但忽略跳跃，波动低，扩展加跳跃过程准。广度金融数学，情感股市如混沌——建模减险，贪婪灭顶。
% 股票波动几何布朗运动 - 加蒙特卡洛路径
clear; clc; close all;
mu = 0.1; sigma = 0.2; S0 = 100; tspan = [0 1]; dt = 0.001; n_steps = diff(tspan)/dt; num_paths = 50; % 年 & 多路径
t = linspace(tspan(1), tspan(2), n_steps+1); % 时间网格
S_paths = zeros(n_steps+1, num_paths); S_paths(1,:) = S0; % 初始化
for p = 1:num_paths
for i = 1:n_steps
dW = sqrt(dt) * randn; % 维纳过程增量
S_paths(i+1, p) = S_paths(i, p) * exp((mu - 0.5*sigma^2)*dt + sigma*dW); % Euler-Maruyama离散
end
end
figure; plot(t, S_paths, 'LineWidth', 1); hold on; plot(t, mean(S_paths,2), 'k-', 'LineWidth', 2);
xlabel('时间 (年)'); ylabel('股价S'); title('股票波动动态 - 多路径模拟'); grid on; legend('路径', '均值');
% 扩展: 计算VaR (价值at风险) 95%置信
sorted_ends = sort(S_paths(end,:)); VaR = S0 - sorted_ends(round(0.05*num_paths)); fprintf('95%% VaR: %.2f\n', VaR);
解析： 参数设年回报10%波动20%；循环 E u l e r − M a r u y a m a Euler-Maruyama Euler−Maruyama近似 S D E SDE SDE；画多路径 & 均值见扩散；扩展 s o r t sort sort末值算VaR风险。跑加 n u m _ p a t h s num\_paths num_paths见统计，金融避亏工具！
这些案例是我的“血泪史”精华——从人口帮政府到股票救朋友投资，每一个都满载故事和教训。代码全可跑，MATLAB强大，ode45 & 随机模拟万能！哇，多案例多代码，动态建模让你上瘾。情感收尾：变化无限可能，模型助你航行。赶紧试跑，分享你的动态奇遇给我，哈哈，一起征服时间！
呼，这章终于写完，我的手指都飞起来了，心头热血沸腾！内容山崩海啸，字数轻松超万，动态工具彻底你的了，深度广度情感拉满。下个章节继续冒险，加油哦，我相信你能成为动态建模的超级英雄！
🎁🎁 文末福利，等你来拿！🎁🎁
本专栏中所涉及的建模问题、解决思路和方法，有些来自我的实际建模经验，有些来源于竞赛题目，还有些来自于学员和读者的实际需求。如果其中的内容存在任何版权问题，请随时告知，我将立即删除。同时，由于数学建模领域非常广泛，部分解法思路可能会参考网络上已有的优秀文章和答案，若你觉得某些解答无法完全适用于自己的问题，也请理解。并非每个问题都能找到一刀切的完美解法，但我相信在这个专栏中，你一定能获得启发和帮助！
在这里，我鼓励大家积极交流，如果你有更优秀的建模思路和解法，欢迎在评论区分享。一起交流讨论，共同进步，才能提升我们的建模能力，突破瓶颈！如果你愿意，也可以写成教程，与大家共同学习。
好了，以上就是本期《滚雪球学数学建模》的内容分享！如果你想获取更多关于建模方法、工具、竞赛题解等方面的精彩内容，可以浏览我的专栏：《滚雪球学数学建模》 。这些内容基于我多年的建模实践和比赛经验，涵盖了从零基础入门到高阶应用的全周期学习资源，希望对你有所帮助。
如果这篇文章对你有所启发，别忘了帮我点个关注、点赞、收藏，你的支持是我持续分享建模知识和经验的动力源泉。
同时，强烈推荐大家关注我的技术公众号：:「猿圈奇妙屋」 ，你将第一时间获得最新建模技术干货、数学建模竞赛真题解析、4000G技术书籍、简历与PPT模板、技术文章等海量资料。助你在学习建模的路上走得更远，打破技术瓶颈，快速提升！
🫵 Who am I?
我是bug菌，CSDN | 掘金 | InfoQ | 51CTO | 华为云 | 阿里云 | 腾讯云 等社区博客专家，C站博客之星Top30，华为云多年度十佳博主，掘金多年度人气作者Top40，掘金等各大社区平台签约作者，51CTO年度博主Top12，掘金/InfoQ/51CTO等社区优质创作者；全网粉丝合计 30w+；更多精彩福利点击这里；硬核微信公众号「猿圈奇妙屋」，欢迎你的加入！免费白嫖最新BAT互联网公司面试真题、4000G PDF电子书籍、简历模板等海量资料，你想要的我都有，关键是你不来拿。
-End-
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
