蹦极数学模型求解

蹦极的数学建模及其龙格-库塔法求解方法 信息学院计算机系 05 级硕士研究生 赵也非 51051201062 商学院情报学 05 级硕士研究生 周自力 51050500186 统计系 05 级硕士研究生 孙晶 51050625004 摘要 本文通过参照中给出的数据，对蹦极者在蹦极过程受到重力，拉力，空气 阻力等受力分析，依据牛顿第二定律，将这种现实生活中连续状态的非线性系统 进行建模，得到一个完整的蹦极数学模型。该模型表现为蹦极者位置 x对下落时 间 t的二阶常微分方程。然后利用 Matlab 编程，采用龙格-库塔法方法，完成了 赛题中所有问题。全文的分析思路如下： 首先，求解空气阻力与速度的关系。题中给出了一组空气阻力和速度的实测 数据，通过程序 BengJi NiHe.m，进行多项式曲线拟合，发现空气阻力和速度符 合二次多项式，求出了二次多项式的系数，验证了该二次多项式具有良好的拟合 效果。 然后，对蹦极者受力分析，发现这是典型的具有连续状态的非线性系统。建 立二维空间坐标模型，并令蹦极者位置为 X.根据牛顿第二定律，列出蹦极模型 的数学表达式，得到蹦极者下落位置 x对下落时间 t的二阶常微分方程。 为简化计算，决定采用计算机对蹦极数学模型进行数值计算和系统仿真。因 为 MATLAB 只能解一阶常微分方程，所以先手工把上面的二阶常微分方程转化成 一阶常微分方程，再采用计算机求解。通过对 Matlab 中不同的龙格-库塔法方法 进行分析后，发现 ode23 方法最适合求解具有连续状态的非线性系统，且精度符 合要求。 因此，程序（BengJi.m, BengJi_Sub.m）中使用 ode23 方法，对蹦极数学模 型进行数值计算和系统仿真。并得出了系统要求的数值解和系统仿真图表。通过 极值和折半法，求出在蹦极绳弹性系数 k=5 时，蹦极者有最大刺激，即在安全的 情况下最接近湖面。此情况下，脚踝受到的最大拉力为 670 磅，蹦极者的最大速 度为 105.1469 英尺/秒，蹦极者反弹回来离起跳点的最短距离为 69.7566 英尺， 并给出了系统仿真图。 将蹦极系统的理论数值解和仿真数值解进行比对验证，误差分析，发现系统 仿真结果符合实际，本数学模型可以客观正确地反映蹦极过程。 最后，论述了此模型的优缺点，讨论了模型的改进，列出了相关参考文献和 术语。 关键字 数学建模 MATLAB 非线性系统 常微分方程 龙格-库塔法方法 多项式曲线拟合 数值计算 系统仿真 蹦极 1 一. 问题重述 蹦极者脚踝上扣有 L = 160 英尺的蹦极绳（蹦极绳的弹性系数 k 未知），空 气阻力与速度有关，给出下列空气阻力和速度的实测数据。 速度: 英尺/秒 10 20 30 40 50 60 70 80 90 100 110 阻力: 磅英尺/秒2 7 24 51 88 135 192 259 336 423 520 627 已知当地的重力加速度g = 32 英尺/秒2 ，脚踝所受的拉力不能超过 800 磅， 若起跳点离湖面 300 英尺，不让蹦极者碰到水面。蹦极者身高H=6 英尺, 质量 m=160 磅, 且蹦极绳的弹性系数是常数。 求：1 蹦极绳的弹性系数 k 为多少时,在保证蹦极者安全的前提下（脚踝受 拉力不超过 800 磅），能让蹦极者得到最大刺激（最接近水面，而又不进入水中）。 记算出来的 K为 K1. 2 K1 下蹦极者的最大速度(英尺/秒) 3 K1 下蹦极者脚受到最大的拉力(磅) 4 K1 下蹦极者反弹回来时离起跳点可能达到的最短距离(英尺) 二. 模型假设 1. 通过给出的一组空气阻力和蹦极者下落速度的实测数据，根据拟合方法， 得出空气阻力和蹦极者下落速度的数关系。 2. 蹦极者在下落过程中，受到重力，空气阻力，蹦极绳拉力的作用，不考 虑其他作用力的影响。 3. 蹦极绳的弹性系数 k是一固定常数，且蹦极绳不会断裂。 4．蹦极者起跳时候无初速度 5. 蹦极者在起跳到头顶碰到湖面之前为一个质点，不考虑高度和翻转等运 动。 6. 蹦极者得到最大刺激指的是在最能接近水面时 三. 符号定义 m 人的质量 g 重力加速度 k 绳的弹性系数 x 蹦极者的位置 . x 位置 x对时间 t的一阶微分，即蹦极者的速度 .. x 位置 x对时间 t的二阶微分，即蹦极者的加速度 F resis ( ) 空气阻力，与速度 有关 . x . x b(x) 绳子的拉力，与蹦极者 x所在位置有关 四. 模型建立 2 1.建立坐标体系 对蹦极场所建立坐标系，以起跳点下 160 英尺为坐标原点，坐标原点以下为 正，坐标原点以上为负，x为蹦极者在此坐标中的位置。考虑到绳是系在脚踝上， 人身高 6英尺，x的有效范围应该在[-160 英尺，134 英尺]以内，才能保证蹦极 者在起跳点起跳，不碰到湖面。 0 -160 140 图 1 X 坐标系 下落过程，在绳子绷紧之前，人向上受到空气阻力，向下受到重力的作用。 在绳子绷紧之后，人向上受到空气阻力和绳子的拉力，向下受到重力的作用 人的脚踝扣有弹性系数为 k 的蹦极绳，依据胡克定律，蹦极绳的拉力 b(x) 和蹦极者位置 x的函数关系为： b(x) = ⎩⎨ ⎧ <= >− 0,0 0, x xkx F resis ( )为空气阻力，与速度有关。将在后面的章节依据给出的 11 组数据 找出它们的关系。 . x 对人进行受力分析可以，可得到下图，mg 为人的重力 b(x) F ( ) resis . x mg 3 图 2 蹦极者受力分析 2 模型分析 在蹦极者下落过程中，蹦极者处于失重状态下。按照牛顿第二定律，质量乘 以加速度等于合外力 m =F （公式 1） .. x 由公式 1以及上面的推导可得出蹦极系统的数学模型为： m = mg - F ( )-b(x) （方程 1） .. x resis . x 这是典型的具有连续状态的非线性系统。 依据假设，蹦极者的起始速度为零， (0) = 0 . x 依据假设，蹦极者在起跳时视为质点，蹦极者的起始位置为蹦极绳的长度， 即 x(0) = -160 五. 模型求解 1．空气阻力与速度之间的关系 F resis ( ) . x 已知给出了一组空气阻力与速度的实测数据，用图形分析得出的结果，如图 3所示： 阻力速度关系 0 100 200 300 400 500 600 700 0 50 100 150 速度 阻 力 阻力: 磅英尺/秒2 图 3 阻力速度关系图 看上去象二次曲线，可以考虑用一个二次多项式进行曲线拟合。 在MATLAB里（BengJi_NiHe.m），对空气阻力和速度进行曲线拟合，得到的结 果图 4所示： 4 图 4 因此得出，空气阻力与速度的关系为： F ( ) = 0.05| | ＋ 0.2 （方程 2） resis . x . x . x . x 2 求解整个蹦极系统的数学模型 将方程 2代入方程 1得到： m = mg + b(x)– 0.05| | - 0.2 b(x) = （方程组 3） .. x . x . x . x ⎩⎨ ⎧ ≤ >− 0,0 0, x xkx 初始条件： x(0) = -160； (0) = 0 . x 注意到 方程组 3 是关于 x的二阶微分方程，需要把上面的二阶微分方程化 为一阶微分方程 ，以方便在 MATLAB 对中求解 令： Y = ⎥⎦ ⎤⎢⎣ ⎡ . x x ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ≤⎥⎦ ⎤⎢⎣ ⎡+⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −−= >⎥⎦ ⎤⎢⎣ ⎡+⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ −−−= 0, 0 ||05.02.00 10 0, 0 ||05.02.0 10 . . . . x g Yx mm Y x g Yx mmm kY 初始条件变为：Y(0) = = ⎥⎦ ⎤⎢⎣ ⎡ )0( )0( . x x ⎥⎦ ⎤⎢⎣ ⎡− 0 160 得到上面的一阶常微分方程后，就可以在 MATLAB 中用龙格-库塔法方法进行 求解，所对应的 MATLAB 命令为 ode（Ordinary Differential Equation 的缩写）。 ode 有 ode23, ode45 等不同方法用于求解不同类型的微分方程。 因为本系统是一个典型的具有连续状态的非线性系统，精度在 10 以内， 所以采用ode23 方法进行求解。程序（ 3− BengJi_Sub.m，BengJi.m） 由已知条件得：对于不同规格的蹦极绳每拉长一英尺拉力增加的范围为 4 至 6磅。蹦极绳的弹性系数 k的取值范围在 4~6 磅/英尺。 湖面的坐标为 140 英尺，所以如果 k选取得当，则蹦极者位置 x的最大值应 该是 134 英尺（因为湖面坐标 140 英尺，且人的身高 6英尺）。 5 首先用两个极值进行试探性求解 当 k = 4 时，如图 4所示，位置 x的最大值在 150 英尺，这时蹦极者落在湖 面以下 10 英尺，需要增大弹性系数 k。 图 4 K=4 的模拟结果 当 k = 6 时，如图 5所示，位置 x的最大值在 115 英尺，这时蹦极者离湖面 还有 40 英尺，没有碰到湖面，需要减小弹性系数 k。 6 图 5 K=6 的模拟结果 弹性系数 K与蹦极者的最大下落距离为单向递减的关系，因此合适的弹性系 数应可以使用折半查找(Binary Search)法（又称二分查找）。 第一次折半，K=(6+4) /2 = 5 时，位置 x的最大值在 134 英尺， 符合实际 情况。所以不用再继续折半查找，因此选取弹性系数 k = 5 磅/英尺。 运行 MATLAB 源程序 BengJi.m，求得蹦极者受到最大拉力，最大速度，离 起跳点的最小距离，得到的数值结果如图 6 所示： 7 图 6 蹦极模型数值解 在蹦极绳的弹性系数 k = 5 时，在 MATLAB 对系统模型进行数值计算得： 1. 他的脚踝受到的最大拉力(磅) ---- F = 670 磅 2. 他的最大速度(英尺/秒) ---- vMax = 105.1469 英尺/秒 3. 他反弹回来时离起跳点可能达到的最短距离(英尺) ---- minL = 69.7566 英尺 六. 验证与分析 当弹性系数k = 5 时，运行MATLAB源程序BengJi.m，系统仿真结果如下图所 示： 图 7 K=5 的模拟结果 其中，横坐标是时间 t，从 0~200 秒。 纵坐标是蹦极者以起跳点下面 160 英尺为坐标原点（因为绳长是 160 英尺）， 坐标原点以上的位置为负值，坐标原点以下为正值。所以蹦极者的起跳点坐标是 -160 英尺，湖面坐标是 140 英尺。 从上图可以看到，蹦极者的起跳位置是-160 英尺，坐标最大位置是 134 英 尺（湖面坐标 140 英尺，人的身高 6 英尺），200 秒之后，蹦极者的坐标位置趋 近于 32 英尺（因为拉力等于重力），这和现实是非常吻合的。 8 七. 模型的优缺点及其改进 1. 优点 1) 求解速度与阻力的关系时，我们拟合到了与测试数据非常符合的模型， 提高了计算结果的精度，这就提高了以后模型最终结果的准确性。 2) 模型的形式是简单，在物理学上下落的过程考虑到受力不同，要分为 4~5 个阶段。而利用向量，微分方程等建立一个形式简单，比较完美的模型。 3) 有很强的操作性和推广性，蹦极者的质量，绳长，蹦极的高度，等都是 可以改变的，改变初值，输入程序即可得到最后不同的结果。 4) 利用 Matlab 编程，使用了龙格-库塔(Runge-Kutta)方法。借助计算机大 大简化了手工计算的工作量 2. 缺点 模型与真实解有一定的差距。现实中，可能还需要考虑风力以及人的翻 转等影响。在此模型中没有予以考虑。如果可以在数据量充足的情况下， 模拟出风力，翻转等影响，模型会更加完美。 八. 参考文献 [1] 刘卫国 陈昭平 张颖，MATLAB 程序设计与应用，北京：高等教育出版社， 2002 年。 [2] 刘兵团，数学建模基础，北京: 清华大学出版社 北京交通大学出版社，2004 年。 [3] 苏晓生，掌握 MATLAB 6.0 及其工程应用，北京：科学出版社 2002 年。 [附录 1] 相关术语 龙格-库塔方法 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法 精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。同前几种算法一样，该算 法也是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： yi＋１＝yi＋h*K1 K1=f(xi,yi) 当用点 xi 处的斜率近似值 K1 与右端点 xi＋１处的斜率 K2 的算术平均值作为平均斜率 K*的近似值，那么就会得到二阶精度的改进欧拉公式： yi＋１＝yi＋h*( K1+ K2)/2 K1=f(xi,yi) K2=f(xi+h,yi+h*K1) 依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值 K1、K2、……Km，并用他 们的加权平均数作为平均斜率 K*的近似值，显然能构造出具有很高精度的高阶计算公式。 9 经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格 －库塔算法： yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 折半查找法： 折半查找法也称为二分查找法，它充分利用了元素间的次序关系，采用分治策略，可在 最坏的情况下用 O(log n)完成搜索任务。它的基本思想是，将 n 个元素分成个数大致相同的 两半，取 a[n/2]与欲查找的 x 作比较，如果 x=a[n/2]则找到 x，算法终止。如果 xa[n/2]， 则我们只要在数组 a 的右半部继续搜索 x。 牛顿第二定律： 物体的加速度跟物体所受的合外力 F 成正比，跟物体的质量成反比，加速度的方向跟合 外力的方向相同。 胡克定律： 在材料的线弹性范围内，固体的单向拉伸变形与所受的外力成正比；也可表述为：在应 力低于比例极限的情况下，固体中的应力σ与应变ε成正比，即σ＝Εε，式中 E 为常数， 称为弹性模量或杨氏模量 10