数学模型与python基础——微分方程数学模型（1）

最新推荐文章于 2024-09-14 08:25:16 发布

一无所知的汪汪

最新推荐文章于 2024-09-14 08:25:16 发布

阅读量1.5k

点赞数 1

分类专栏：数学建模文章标签： python 数学建模算法数学

本文链接：https://blog.csdn.net/qq_34345989/article/details/121903001

版权

数学建模专栏收录该内容

3 篇文章 0 订阅

订阅专栏

本文通过数学模型分析火箭发射过程，解释为何必须使用多级火箭来将卫星送入高空轨道。首先，建立了单级火箭模型，利用微分方程和Python的sympy库求解，得出火箭速度与火箭质量比的关系，表明一级火箭无法满足要求。然后，构建理想火箭模型，同样利用微分方程求解，揭示多级火箭如何通过逐级抛弃无用结构来提高速度。最后，通过实例计算展示三级火箭相较于二级火箭的优势，从而证明多级火箭在技术上的必要性。

摘要由CSDN通过智能技术生成

数学模型与python基础（3.1）

第一章讲了什么是数学模型，介绍了传染病学的sir模型。那么我们就趁热打铁从微分方程相关的数学模型以及相关的数学计算过程开始介绍。

微分方程数学模型

因为有很多系统虽然很难直接描述变量之间的关系，但是很容易描述变量和其微小增量或变化率的关系式，因此可以通过微分方程构件系统的数学模型。建模时，我们需要根据建模目的和问题分析做出简化假设，并且按照内在规律或类比法建立微分方程。

数学建模是有固定步骤的，必须要按照相应的步骤来才能高效地解决问题。步骤如下：
（1）问题的提出

（2）模型假设

（3）模型建立与求解

（4）结果分析

（5）模型评价与推广

数学建模的论文也应该以此作为行文顺序撰写，以下将以几个实际的例子来说明，不过由于具体的结果分析和模型评价与推广要基于具体问题和大量图表，在接下来的例子中就不涉及了。

火箭发射

（1）问题的提出

采用运载火箭把人造卫星发射到高空轨道上运行，为什么不能用一级火箭而必须用多级火箭系统？

（2）模型假设

A.卫星轨道为过地球中心某一平面的圆,卫星在此轨道上以地球引力做为向心力绕地球做平面圆周运动

B.地球是固定于空间中的均匀球体，其它星球对卫星引力忽略不记。地球半径为R，中心为O，地球质量看成集中于球心(认为地球是均匀球体)。

曲线 $C$ 为地球表面， $C^{'}$ 为卫星轨道，其半径为 $r$ ，卫星质量为 $m$ 。

根据牛顿定律,地球对卫星的引力为:
$F=G\times \frac{m}{r^2}$
其中G为引力常数，可根据卫星在地面的重量算出，
$\frac{Gm}{R^2}=mg,G=gR^2$
因此
$F=mg\times (\frac{R}{r})^2$
由于卫星所受引力为其做匀速圆周运动时的向心力
$F=m\times \frac{v^2}{r}$
卫星速度为
$v=R\times \sqrt{\frac{g}{r}}$
假定需要将卫星送入600km的高空轨道，即 $r - R = 600 k m$ ，取地球半径 $R = 6400 k m$ ，取重力加速度 $g=9.81m/s^2$ ，那么计算可得 $v = 7.6 k m / s$

C.火箭在喷气推动下做近似直线运动，空气阻力不记

D.火箭喷出同种气体相对火箭的速度是不变的，假设目前的技术条件，火箭喷出气体相对火箭的速度 $u = 3 k m / s$

（3）模型建立与求解

单级火箭模型

由于燃料燃烧会喷出气体，所以火箭质量不是始终恒定不变的，假定火箭喷出气体相对火箭的速度为 $u$ ，则气体相对地球运动的速度为 $v (t) - u$

对火箭飞行过程中的质量变化进行泰勒展开可得
$m(t+\Delta t)-m(t)=\frac{dm}{dt}\Delta t+o(\Delta t^2)$
$v(t+\Delta t)-v(t)=\frac{dv}{dt}\Delta t+o(\Delta t^2)$

由 $\Delta t$ 前后动量守恒可得
$v(t)=m(t+\Delta t) v(t+\Delta t)-\left(\frac{d m}{d t} \Delta t+o\left(\Delta t^{2}\right)\right)(v(t)-u)$

带入泰勒展开并忽略高阶小项
$v(t)=[m(t)+\frac{dm}{dt}\Delta t][v(t)+\frac{dv}{dt}\Delta t]-\frac{d m}{d t} \Delta tv(t)+\frac{d m}{d t} \Delta tu$

$v(t)+m(t)\frac{dv}{dt}\Delta t+\frac{dm}{dt}\Delta tv(t)+\frac{dm}{dt}\Delta t\frac{dv}{dt}\Delta t-\frac{d m}{d t} \Delta tv(t)+\frac{d m}{d t} \Delta tu$

忽略 $\Delta t^2$ 项
$0=m(t)\frac{dv}{dt}\Delta t+\frac{dm}{dt}\Delta tv(t)-\frac{d m}{d t} \Delta tv(t)+\frac{d m}{d t} \Delta tu$

$\frac{d v}{d t}\Delta t =-u \frac{d m}{d t}\Delta t$

$\frac{d v}{d t} =-u \frac{d m}{d t}$

解微分方程。我知道有很多人已经忘掉，或者也没怎么会过微分方程的解法，这没关系，我们有python。

python有一个sympy科学计算库，具有大量符号计算的功能，解方程、求导、求极限、表达式求值之类的都能干。别的功能用到再说，这篇讲一下它怎么解常微分方程。

求解这个微分方程的代码如下

import sympy as sy
u, t, m0, v0 = sy.symbols('u, t, m_0, v_0')  # 为变量定义符号
v = sy.Function('v')  # 定义函数
m = sy.Function('m')  # 定义函数
# 定义微分方程,其中diff(v(t),t,1)表示v(t)对t求1阶导
ode = m(t)*sy.diff(v(t), t, 1)+u*sy.diff(m(t), t, 1)
# 上式中的微分方程表示m(t)*dv(t)/dt+u*dm(t)/d(t)=0
ode_sol = sy.dsolve(ode, v(t))  # 用dsolve函数解微分方程,且求解的是v(t)的表达式
sy.pprint(ode_sol)  # pretty print，美化打印，输出v(t)的通解
ics = {v(0): v0, m(0): m0}  # 定义微分方程的初始条件
# 将微分方程的初始条件带入通解，构建关于常数的方程，并用solve函数求解
# subs为替换函数用值或表达式替换原式中的值或表达式
C_sol = sy.solve(ode_sol.subs(t, 0).subs(ics))
sy.pprint(ode_sol.subs(C_sol[0]))  # 输出v(t)的特解

输出如下

v(t) = C₁ - u⋅log(m(t))
v(t) = u⋅log(m₀) - u⋅log(m(t)) + v₀

第一个输出为v(t)的通解，第二个输出为v(t)的特解，所以我们知道在该数学模型中
$v(t)=v_{0}+u \cdot \ln \left(\frac{m_{0}}{m(t)}\right)$
通过该式我们得知，当 $v_0$ 和 $m_0$ 一定的情况下，火箭速度 $v (t)$ 由火箭喷出气体相对火箭运动速度 $u$ 以及初始质量与当前质量的比值 $\frac{m_0}{m(t)}$ 。因此提升火箭速度可以通过提升u或者减少m(t)达到，由于 $u$ 不是本题所讨论的对象，因此仅对火箭质量m进行建模。

火箭质量 $m (t)$ 由三部分构成，即火箭负载 $m_p$ ，燃料质量 $m_F$ ，结构质量 $m_s$

假如用一级火箭发射卫星，当燃料耗尽时，即 $m_F=0$ 时， $m(t)=m_p+m_s$ ，火箭初始速度 $v_0=0$ ，则
$\cdot \ln \left(\frac{m_{0}}{m_{p}+m_{s}}\right)$
为了满足在 $u = 3 k m / s$ 的情况下让 $v > 7.6 k m / s$ ，需要 $\ln \left(\frac{m_{0}}{m_{p}+m_{s}}\right) \ge 2.53$ 了解技术条件发现， $\frac{m_{s}}{m_{F}+m_{s}} \geq \frac{1}{9}$

假设
$m_{s}=\lambda\left(m_{F}+m_{s}\right) = \lambda\left(m_{0}-m_{F}\right)$
则
$\cdot \ln \frac{m_{0}}{\lambda m_{0}+(1-\lambda) m_{p}}$
假设火箭无负载，取 $\lambda = \frac{1}{9}$ ，则 $\cdot ln\frac{1}{\lambda} = u \cdot ln9 \approx 6.6km/s<7.6km/s$ 无法送入距地600km的高空轨道

由此我们可以知道，仅用一级火箭在空载的情况下都无法进入高空轨道，更无法将卫星送入了，但是建模到这一步还并没有完成，还需要论证为什么多级可以。

理想火箭模型

多级火箭的模型相对而言复杂，因此先构件理想火箭模型，认为火箭每时每刻都会抛弃一部分无用结构来减轻自身重量。假设每时每刻理想火箭损失的质量中有一部分为无用的结构质量，还有一部分为消耗的燃料质量，这两部分的比值为 $\frac{\alpha}{1-\alpha}$ ，因此根据动量定理可得下式
$v(t)=m(t+\Delta t) v(t+\Delta t)-\alpha \frac{d m}{d t} \Delta t \cdot v(t) -(1-\alpha) \frac{d m}{d t} \Delta t \cdot(v(t)-u)+o(\Delta t)$
带入泰勒展开并忽略高阶小项
$0=m(t)\frac{dv}{dt}\Delta t+\frac{dm}{dt}\Delta tv(t)+\frac{dm}{dt}\Delta t\frac{dv}{dt}\Delta t-\alpha \frac{d m}{d t} \Delta t \cdot v(t) -(1-\alpha) \frac{d m}{d t} \Delta t \cdot(v(t)-u)$
消去 $\Delta t$ 且忽略 $\Delta t^2$ 项
$0=m(t)\frac{dv}{dt}+\frac{dm}{dt}v(t)-\alpha \frac{d m}{d t} \cdot v(t) -(1-\alpha) \frac{d m}{d t} \cdot(v(t)-u)$

$-m(t)\frac{dv}{dt}=\frac{dm}{dt}v(t)-\alpha \frac{d m}{d t} \cdot v(t) -\frac{d m}{d t} \cdot(v(t)-u)+\alpha\frac{d m}{d t} \cdot(v(t)-u)$

$-m(t)\frac{dv}{dt}= \frac{d m}{d t} \cdot u-\alpha\frac{d m}{d t} \cdot u$

即
$\frac{d v(t)}{d t}=(\alpha - 1) \frac{d m}{d t} \cdot u$
修改刚才的代码，求解

import sympy as sy
u, t, m0, v0, a = sy.symbols('u, t, m_0, v_0, a')  # 为变量定义符号
v = sy.Function('v')  # 定义函数
m = sy.Function('m')  # 定义函数
# 定义微分方程,其中diff(v(t),t,1)表示v(t)对t求1阶导
ode = m(t)*sy.diff(v(t), t, 1)+(1-a)*u*sy.diff(m(t), t, 1)
# 上式中的微分方程表示m(t)*dv(t)/dt+(1-a)*u*dm(t)/d(t)=0
ode_sol = sy.dsolve(ode, v(t))  # 用dsolve函数解微分方程,且求解的是v(t)的表达式
sy.pprint(ode_sol)  # pretty print，美化打印，输出v(t)的通解
ics = {v(0): v0, m(0): m0}  # 定义微分方程的初始条件
# 将微分方程的初始条件带入通解，构建关于常数的方程，并用solve函数求解
# subs为替换函数用值或表达式替换原式中的值或表达式
C_sol = sy.solve(ode_sol.subs(t, 0).subs(ics))
sy.pprint(ode_sol.subs(C_sol[0]))  # 输出v(t)的特解

输出为

v(t) = C₁ + u⋅(a - 1)⋅log(m(t))
v(t) = -a⋅u⋅log(m₀) + u⋅(a - 1)⋅log(m(t)) + u⋅log(m₀) + v₀

即
$v(t)=v_0 + (1-\alpha) u \ln \frac{m_{0}}{m(t)}$
由 $v (0) = 0$
$v(t)=(1-\alpha) u \ln \frac{m_{0}}{m(t)}$
在这种模型下，想要达到 $v (t)$ 只需要一个足够大的 $m_0$ 即可

多级火箭模型

理想火箭模型确实是很好的假设，但是现在的技术无法实现，因此构件多级火箭模型。假定多级火箭的工作模式为火箭由每级质量为 $m_i$ 的多级火箭构成，当第i级燃料烧尽时，第i +1级火箭立即自动点火，并抛弃已经无用的第i级。假设每级火箭的总质量比结构质量的比值 $\lambda$ 是相同的，则每级火箭的结构质量为 $\lambda m_i$ ，燃料质量为 $(1-\lambda)m_i$ 认为火箭喷出燃料的速度始终相对火箭速度为u，且假设每一级的初始质量与其负载质量（该级火箭的负载包括所有火箭的负载，如卫星等，以及未工作的火箭质量）之比为 $k_i$ 。若多级火箭可行的话，求解k使得 $\frac{m_0}{m_p}$ 最小，则是我们的多级火箭方案。

定义
$\begin{aligned} w_{1}=& m_{0}=m_{1}+m_{2}+\cdots+m_{n}+m_{p} \\ w_{2}=& m_{2}+\cdots+m_{n}+m_{p} \\ & \cdots \\ w_{n+1} &=m_{p} \end{aligned}$
则根据上文定义
$k_{1}=\frac{w_{1}}{w_{2}}, \cdots, k_{n}=\frac{w_{n}}{w_{n+1}}$
第i级的末速度为
$\begin{array}{c} v_{1}=u \ln \frac{m_{1}+m_{2}+m_{p}}{\lambda m_{1}+m_{2}+m_{p}}=u \ln \frac{k_1+1}{\lambda k_1+1} \\ v_{2}=v_{1}+u \ln \frac{m_{2}+m_{p}}{\lambda m_{2}+m_{p}}=u \ln \frac{k_1+1}{\lambda k_1+1}+ u \ln \frac{k_2+1}{\lambda k_2+1}=uln(\frac{k_1+1}{\lambda k_1+1}\cdot \frac{k_2+1}{\lambda k_2+1})\\ \dots\\ v_{\text {末 }}=u \ln \left(\frac{w_{1}}{\lambda m_{1}+w_{2}} \cdots \frac{w_{n}}{\lambda m_{n}+w_{n+1}}\right) \end{array}$

将速度和 $\frac{m_0}{m_p}$ 均写为k的表达式
$\begin{array}{c} m_{1}=w_{1}-w_{2}, \quad m_{2}=w_{2}-w_{3}, \cdots, \quad m_{n}=w_{n}-w_{n+1} \\ v_{\text {末 }}=u \ln \left(\frac{k_{1}}{\lambda\left(k_{1}-1\right)+1} \cdots \frac{k_{n}}{\lambda\left(k_{n}-1\right)+1}\right) \\ \frac{m_{0}}{m_{p}}=k_{1} k_{2} \cdots k_{n} \end{array}$
则问题的数学表述为，求解最小的 $k_1k_2...k_n$ 使得
$\frac{k_{1}}{\lambda\left(k_{1}-1\right)+1} \cdots \frac{k_{n}}{\lambda\left(k_{n}-1\right)+1}=e^{\frac{v_{\text {末 }}}{u }}$
推知当各级k均相同时， $\frac{m_{0}}{m_{p}}=k_{1} k_{2} \cdots k_{n}$ 最小

假如是二级火箭的话，即
$\frac{m_1}{m_2+m_p}=\frac{m_2}{m_p}=k$
第一级、第二级火箭燃烧完时，速度分别为
$\begin{array}{c} v_{1}=u \ln \frac{m_{1}+m_{2}+m_{p}}{\lambda m_{1}+m_{2}+m_{p}}=u \ln \frac{k+1}{\lambda k+1} \\ v_{2}=v_{1}+u \ln \frac{m_{2}+m_{p}}{\lambda m_{2}+m_{p}}=2 u \ln \frac{k+1}{\lambda k+1} \end{array}$
若 $\lambda = 0.1$ ， $u = 3 k m / s$ 且综合考虑空气阻力等因素之后，认为我们需要让第二级火箭速度 $v_2=10.5km/s$
$\begin{array}{c} 6 \ln \frac{k+1}{0.1 k+1}=10.5 \\ k=11.2, \\ \frac{m_{0}}{m_{p}}=\frac{m_{1}+m_{2}+m_{p}}{m_{p}}=(k+1)^{2} \approx 149 \end{array}$
若为3级火箭，则
$v_{3}=3 u \ln \frac{k+1}{\lambda k+1}$
同样的条件下求得 $k = 3.25$ ， $\frac{m_0}{m_p}=77$