image.png
工程与科学数值方法的matlab实现（第4版）
[美]Steven C.Chapra 著
林赐 译
长久以来，人们都习惯用解析法研究问题。只有解析法无效的时候，才考虑用数值分析方法来解决问题。所以它一直被当做备选项，不被重视。
数值分析方法是一门实践性很强的学科，但是其计算繁琐、冗长，靠人工计算计算量大让人望而生畏。随着计算机的发展，尤其出现了一批优秀的数值分析软件比如MATLAB软件，完美的解决了数值分析的缺点，使数值分析方法被广泛使用。成为与理论分析、科学实验并列的第三种科学研究手段。
image.png
第1部分建模、计算机与误差分析
数值方法（numerical method）使用公式表示数学问题以便可以用算数和逻辑运算解决这些问题的技术。数字计算机擅长于执行算术和逻辑这类运算。所以数值方法又称为计算机数学（computer mathematics）。
数值方法学习的理由：
（1）数值方法能够极大地覆盖所能解决的问题类型，增强问题求解技巧；
（2）数值方法可以让用户更加智慧地使用“封装过的”软件；
（3）很多问题不能用封装的程序解决，可以自我编程解决问题；
（4）数值方法是学习使用计算机的有效载体；
（5）数值方法提供了一个能够增强对数学理解的平台。
第1章 数学建模、数值方法与问题求解
本章目标：
（1）学习如何基于科学原理建立数学模型，进而对简单物理系统的行为进行仿真；
（2）理解数值方法如何提供一种方式以便在数字计算机上求得问题的解；
（3）对于工程学科中使用的各种模型，理解背后的守恒律，正确评价这些模型稳态解和动态解之间的差异；
（4）学习本书中设计的不同类型的数值方法。
提出问题：预测蹦极运动员自由落体阶段的速度。
image.png
解：根据牛顿第二定律就可以建立速度的关系式：
image.png
image.png
根据空气流体力学有：
image.png
整理后有最终的表达式为：
image.png
其中，v：垂直向下速度(m/s)
t：时间(s)
g：重力加速度(≈9.81m/s2)
Cd：集总阻力系数(kg/m)
（集总是指它的大小还取决于其它一些因素，比如运动员的表面积和流体密度）
m：为蹦极运动员的体重(kg)
1.1 一个简单的数学模型
数学模型（mathematical model）可以广义地定义为表达物理系统或物理过程本质特性的公式或方程。
广义数学模型可表示为如下函数关系式：应变量=f（自变量，参数，强制函数）
其中，应变量（dependent variable）是物理系统的特征，反映了系统的行为或状态；
自变量（independent variable）通常为维度，如时间和空间，通过维度确定系统的行为；
参数（parameter）是系统特性或构成的反映；
强制函数（forcing function）是作用在系统上的外部影响。
（1.1）式的特解为：（具体求解过程见习题1.1）
image.png
其中tanh为双曲正切函数。
image.png
例1.1蹦极问题的解析解：
设运动员m=68.1kg，空气阻力Cd=0.25kg/m，g=98.1m/s2，计算t从0~12s内自由落体过程中运动员的速度。
解：将相关数据代入（1.2）式简化后如下式：
V(t)=51.6938tanh(0.18977t)
把时间t=0~12s的值带入上式，列成表1如下：
image.png
继续增加时长到速度不变，如表2：
image.png
image.png
从表2可以看出，经过足够长时间后，下降的速度达到恒定值，即重力与大气阻力处于平衡，此时的速度即为极限速度。根据上表可以看出，极限速度为51.6938.绘制的图形如下：
image.png
式（1.2）又称为解析解或闭型解。因为它与初始微分方程是精确对应的。然而大部分的数学模型还不能精确求解，只能尽力逼近精确解的数值解。
如果（1.1）式无法求得精确解，可以通过下面方法到逼近值：
image.png
上式称为时刻一阶导数的有限差分逼近。代入式（1.1）可得：
image.png
整理后有：
image.png
将（1.1）式带入上式，可以将上式简化为：
image.png
（1.4）式表示新值可以用初始值和斜率及步长有关，如果能确定这三个值就可以确定下一个值。即：新值=上一时刻的值+斜率*步长
这种方法称为欧拉法。
例1.2 蹦极问题的数值解
利用例1.1的数据，用欧拉法计算速度。
image.png
image.png
image.png
image.png
当步长越小的时候，精确解和近似解越逼近。
1.2 工程与科学中的守恒律
守恒律：变化量=增加量-减少量
如果变化量=0，增加量=减少量，即为稳态计算。守恒律也可以应用到化学工程、土木工程、机械工程、电子工程等计算中。
image.png
1.4 案例研究 在以前的学习中，假定阻力为取决于速度平方 ，瑞丽爵士给出了更详细的表达式，可以写为：
image.png
其中，
Fd=阻力(N),
[图片上传失败...(image-2ee934-1585990641587)] =流体密度(kg/m3),
A=垂直于运动方向平面上物体的正面面积(m3),
Cd=无量纲阻力系数，
v=速度方向的单位向量。
假定了湍流条件（即高雷洛数）的关系，集总阻力系数可以表示为：
image.png
即集总阻力系数取决于物体的面积、流体密度和无量纲阻力系数。
注意：
对于速度比较低的情况，物体周围的流动状态呈现出层流，并且阻力与速度之间的关系变成了线性关系。即所谓的斯托克斯阻力。
由于速度平方后，符号和方向丢失。所以修改模型，使其同时适用于向下和向上的速度。
image.png
image.png
待求的微分方程为：
image.png
由此也可以求出位置x(m)：
image.png
利用欧拉法可以将距离方程改成：
image.png
总结： 修正前的速度：
image.png
修正后的速度：
image.png
带入已知条件，v(0)=-40,x(0)=0，t=2s，分别求修正前和修正后的速度和位移数据：
image.png
image.png
从这个图形可以看出，速度修正值和原来值中间有差值。
