数学模型中的反问题逆问题

数学模型中的反问题
向下运动
向上运动
风筝
数学模型竟赛中有很多涉及反问题。

如2010国赛中A题和2011年美赛中A题都涉及反问题。

顾名思义，反问题是相对于正问题而言的。

正问题的定义为：按着自然顺序来研究事物的演化过程或分布形态，起着由因推果的作用。

自然顺序的定义为：不受任何限制和约定俗成的顺序，一般地都认为他们是自然而然的，无须多加解释的。

在一般地语境下，认为这些顺序都是是前提条件的。

如时间顺序、空间顺序、因果顺序，等等。

纯粹的自然顺序的例子是第一，第二，第三这种升序；或者反过来的倒序；约定俗成的例子是上北下南左西右东。

反问题的定义为：根据事物的演化结果，由可观测的现象来探求事物的内部规律或所受的外部影响，由表及里，索隐探秘，起着倒果求因的作用。

可以看出，正、反两方面都是科学研究的重要内容。

但相对正问题，反问题求解难大，计算量大。

许多人知道求解问题的思路，但由于选用计算方法不适当，在几天内求不出计算结果，失去获奖机会。

尽管一些经典反问题的研究可以追溯很早，反问题这一学科的兴起却是近几十年来的事情。

在科学研究中经常要通过间接观测来探求位于不可达、不可触之处的物质的变化规律；生产中经常要根据特定的功能对产品进行设计，或按照某种目的对流程进行控制。

这些都可以提出为某种形式的反问题。

可见，反问题的产生是科学研究不断深化和工程技术迅猛发展的结果，而计算技术的革命又为它提供了重要的物质基础。

现在，反问题的研究已经遍及现代化生产、生活、研究的各个领域。

简单的概括不足以说明问题，我们下面具体介绍一些常见的反问题类型，希望大家能够对它有一个概括的了解.
第一节反问题的例子
例1 物体下落距离L与时间T，正问题是：已知物体的高度，测量下落时间，即t=t(x). 反问题是：已知物体下落时间，求物体的高度,即x=x(t)。

当人们不知道自由落体运动规律x=0.5gT2之前，能用时钟测量物体下落时间，但反过来，给定下落时间，测量物体高度比较难。

对于没有读中学的人，能完成时钟测量物体下落时间的试验。

但给他物体下落时间，测量物体的下落高度是不容易的事情。

例2 年龄与身高。

正问题是，根据年龄T ，每周岁测身高H ，得到身高H 与年龄T 的关系H=H(T). 反问题是：已知身高H ，求年龄T ，即求关系T=H(T).
例3速度V 与轨道形状z=f(x)，其摩擦系数为μ，z 为高度，初始速度为V0，末速度为Ve=V(y=H). 正问题是，已知轨道形状z=f(x)，求末速度为Ve.反问题是：给定末速度为Ve.求轨道形状z=f(x)。

对大学生，正问题能求出来，但反问题有些难。

例4 热传导问题（2013美赛A 题）
设点P 到边界的距离为x, 传热系数为a , 温度T=T(a,x). 正问题是：已知距离x, 传热系数a 的数值, 求温度T 。

如用一维传热公式：
)sin()exp(2ππx t a T -=
反问题是：已知温度T,x ，求传热系数为a ,
例5 光电板问题(2012A 题)
设屋顶面积为D ，光电板长为L ，宽为H 。

正问题是：已知D ，L ，H ，求在屋顶上铺设光电板最大数量N 。

反问题是：已知光电板总铺设面积D*，光电板长L ，宽H ，求屋顶面积。

由上面几个例子，可以在数学上定义正问题为y=f(x)，定义域为D，值域为V。

反问题为x=g(y). 由高等数学可知，若函数f(x)在D上是单调的，则反函数g(y)存在且唯一。

相对正问题而言，反问题计算量大，选用适当的计算方法是成功求解反问题的关键。

因而要求在求反问题之前，要求学生掌握基本的计算方法。

第二节计算方法
2.1方程求根
在数学建模中，求解方程的根是经常遇到的。

常用求根方法有迭代法，二分法，牛顿法，极小值法，一维寻查法，格子法。

2.1.1迭代法
设函数f(x)=x-g(x)有一根x*, 则f(x*)=0, 或x*-g(x*)=0; 或x*=g(x*); 定义求根的迭代公式为：
定理: 若导数g ’的绝对值小于1, 即|g ’|≤L<1, 则迭代收敛。

证：由于x*=g(x*)，则
x k+1-x*=g(x k )-g(x*)=g ’(ξ) (x k -x*) 有
| x k+1-x*|<L| x k -x*|<L 2| x k-1-x*|<L k+1
| x 0-x*|
因为L<1, 则极限L k -->0, 故x k →x*. 证毕。

例 求f(x)=x-x*x 的零点。

解：这里g(x)=x*x, g ’(x)=2x, 则当|x|<0.5时，|g ’(x)|<1, 即|x|<0.5时，迭代公式
. x k+1=x k 2 收敛。

取x0=0.1, 计算得
X 1=x 02=0.12=10
-2 X 2=x 12=(10-2)2=10
-4
…….. 最后求得x k →x*=0. 实际上，我们知道x=0为x=x*x 的解，但它还有一解x=1; 由于|2x|=|2*1|=2, 则用上面迭代公式x=g(x)=x*x 求不出解x=1. 它需要构造另一种迭代公式 ,...2,1,0)(1==+k x g x k k
. x k+1=g(x k)=√x k
容易验证当x=1时，|g’|<1. 取x0=2, 计算得
X1=x00.5=20.5≈1.414
X2=x10.5=(20.5) 0.5=20.25≈1.1189
X3=x30.5=(20.25) 0.5=20.125≈1.090
……..
最后求得x k→x*=1.
由上面例子可知，对同一函数f(x)，它的不同零点对应的迭代公式不同。

2.1.2二分法
在高等数学里，我们已学习下面定理。

定理：设f(a)f(b)<0, f(x)在区间[a,b]上连续可导，则至少有一个(a,b)中的点x*，使
f(x*)=0.
取a0=a, b0=b, x1=(a0+b0)/2, 则点x*属于子区间[a0,x1]，或子区间[x1,b0]。

若属于子区间[a0,x1]，取a1=a0,b1=x1. 否则属于子区间[x1,b0]，取a1=x1,b1=b0. 得到点x*属于子区间[a1,b1]，且b1-a1=0.5(b0-a0)，即区间长度只有原始区间的一半。

类似上面方法，取x2=(a1+b1)/2, 则点x*属于子区间[a1,x2]，或子区间[x2,b1]。

若属于子区间
[a1,x2]，取a2=a1,b2=x2. 否则属于子区间[x2,b1]，取a2=x2,b2=b1. 得到点x*属于子区间[a2,b2]，且b2-a2=0.5(b1-a1)=0.25(b0-a0), 即区间长度只有原始区间的四分之一。

这种方法一直分二去，得点x*属于子区间[ai,bi], 和数列{xi}, 且bi-ai=0.5i(b0-a0)→0, xi→x*.
例. 求f(x)=1-x2在区间[0.5, 2]上的零点。

解这里a0=0.5, b0=2; 有f(a0)=f(0.5)=1-0.52=0.75, f(b0)=f(2)=1-22=-3, 有
f(a0)f(b0)=0.75*(-3)<0, 故在[0.5,2]上f(x)有一零点x*. 取
x1=(a0+b0)/2=(0.5+2)/2=1.25, 有f(x1)=f(1.25)=-0.5625,
f(a0)f(x1)=0.75*(-0.5625)<0, 则在区间[a0,x1]中有零点x*, 故取a1=a0=0.5,
b1=x1=1.25, x2=(a1+b1)/2=0.875, 计算得f(x2)=f(0.875)=1-(0.875)2=0.2343,
F(a1)f(x2)=0.75*0.2343>0, 则零点x*在区间[x2,b1]=[0.875,1.25]中，故取a2=0.875, b2=1.25. 如此计算下去，当bi-ai<ε=0.3时，求得xi=x3=1.0625。

它为x*的近似值。

二分法的计算步骤为：
1)输入a0,b0, 误差限ε；
2)若f(a0)f(b0)>0, 无根，停止计算。

否则转下一步；
3)取x1=(a0+b0)/2, 若f(a0)f(x1)<0, 取a1=a0,b1=x1; 否则取a1=x1, b1=b0;
4)若b1-a1<ε, 输出近似根x*=(a1+b1)/2; 否则a1→a0, b1→b0, 转第三步。

二分法能用图形来说明，其示意图见图2.1, 图中给出了点a0,b0,x1,x2,x3, 它们根据二分法计算。

由图可知，当二分次数增加时，中间点xi相互靠近，收敛于零点x*.
图2.1 二分法示意图
2.1.3 极小值法
定型：若x*为f(x)的零点，则它为F(x)=f2(x)的极小值点。

证：由于F(x)非负，F(x*)=f2(x*)=00=0, 则x*为F(x)的一个极小值点。

我们容易得：
定理：若F(x) =f2(x),F(x*)=0, 则f(x*)=0.
可见，f(x)的零点计算问题能化为极小值计算问题。

它常用一维寻查法求解。

一维寻查法比较简单，它的计算步骤为
1)输入初始点d0, 步长h, 误差ε；
2)计算函数值F(d0),F(d0-h),F(d0+h);
3)若F(d0)<min{F(d0-h),F(d0+h)}; 转第6步。

4)若F(d0)>F(d0-h), 取d1=d0-h; 否则取d1=d0+h;
5)令d1→d0, h→2h, 转第2步。

6)取a=d0-h, b=d0+h, 用二分法求极值点。

二分法求极值点的原理与求根原理类似。

由下面定理给出：
定理：设F(x)在[a,b]上连续，且0≤F(x),若c为[a,b]中的点，且 F(c)<min{F(a),F(b)}, 则F(x)在[a,b]上存在极小值点x*.
上面定理用反证法容易证明。

二分法求极值点的步骤为，取a0=a,b0=b,c0=c=(a0+b0)/2; h=(b-a)/2; 中点x0=(a0+c0)/2; y0=(c0+b0)/2; 若
F(z0)=min{F(a0),F(b0),F(x0),F(c0),F(y0)}, z0为{a0,c0,b0,x0,y0} 中的一点，取
h1=h/2, a1=z1-h1, c1=z1, b1=z1+h1; 容易计算出(b1-a1)=0.5(b0-a0); 即长度只有原区间的一半。

用类似方法，取中点x1=(a1+c1)/2; y1=(c1+b1)/2; 若
F(z1)=min{F(a1),F(b1),F(x1),F(c1),F(y1)}, z1为{a1,c1,b1,x1,y1} 中的点，取
h2=h/4, a2=z2-h2, c2=z2, b2=z2+h2; 有(b1-a1)=0.5(b1-a1)= 0.25(b0-a0); 则长度只有初始区间长度的四分之一。

如此下去，我们得到点ai, ci, bi, 且
(bi-ai)=0.5i(b0-a0)→0. 可以证明，xi→x*为极小值点。

由上面讨论可知，求极小值点分为两步，先求极点所在的区间[a,b], 然后用二分法逐步缩小区间，求出极小值点。

其计算过程可以用图2.2说明。

图中函数F只有一个极小值点。

给定初值d0和步长h, 求出d0+h为最小值，取2h, 计算得d0+2h也为最小，再取4h, 计算得d0+2h也为最.
图2.2 极小值示意图
例. 用极小值法求函数f(x)=1-x*x的零点，x0=1.4, h=0.1.
解. 令F(x)= f2(x*)=(1-x2)2
先用一维寻查法求含有根的区间[a,b]. 计算F(x0-h)=F(1.3)=0.4761;
F(x0)=F(1.4)=0.9216; F(x0+h)=F(1.5)=1.5612; 比较3个数值，x1=1.3时F=0.4761最小。

将步长放大2倍，取h=0.2, 计算F(x1-h)=F(1.1)=0.0441; F(x1)=F(1.3)=0.4761;
F(x1+h)=F(1.5)=1.5612; 比较3个数值，x2=1.1时F=0.0441最小。

再将步长放大2倍，取h=0.4, 计算F(x2-h)=F(0.7)=0.216; F(x2)=F(1.1)=0.0441; F(x2+h)=F(1.5)=1.5612; 比较3个数值，x3=1.1时F=0.0441最小。

因而取a=0.7, b=1.5.
再用二分法求极值点。

取a0=0.7,b0=1.2,c0=0.95, h=0.25; 中点x0=(a0+c0)/2=0.825; y0=(c0+b0)/2=1.075; 计算得F(0.7)=0.216; F(0.825)=0.102; F(0.95)=0.0095;
F(1.075)=0.02421; F(1.2)=0.1936; 当z1=0.95时，函数F(0.95)=0.0095最小。

则取
a1=0.825, c1=0.95, b1=1.075; x1=0.8875; y1=1.0125; 计算得F(0.825)=0.102;
F(0.8875)=0.04509, F(0.95)=0.0095; F(1.0125)=6.328e-4, F(1.075)=0.02421; 给定误差ε=0.1时，若(bi-a1)/4<ε, 输出近似根x*=1.0125.
2.1.4格子法
对于高维问题，格子法是求极值点的常用方法。

它的思想与二分法类似，基本原理为，给定非负的高维函数y=F(X), 初始点X0, 步长h, 将每个坐标分量加上h和减去h, 求最小值y0=min{F(X0),F(X-h),F(X+h)}, 和对应的坐标点X1, 若X1=X0，取h→h/2, 步长减半，否则取h→2h, 步长加倍，再将X1的每个坐标分量加上h和减去h, 求最小值点X2，如此下去，直到步长h<ε为止。

最后Xi为近似最小值点。

例.求方程式组的极小解：x-y(x+y)-1=0; y(x+y)-y-2=0;
解：令
F(x,y)=[x-y(x+y)-1]2+[y(x+y)-y-2]2
取初值点X0=(0,0), 步长h=0.8; 计算得：
F(-0.8,0)=7.24; F(0.8,0)=4.04 F(0,0)=5.0; F(0,-0.8)=3.00; F(0,0.8)=7.355;
可知X1=(0,-0.8)为最小值点，取加倍步长h=1.6, 计算中心得：
F(-1.6,0.8)=20.9; F(1.6,0.8)=4.923; F(0,0.8)=3.0; F(0,-2.4)=83.6; F(0,0.8)=7.32;
则X2=X1=(0,-0.8)为最小值点，取减半步长h=0.8,继续计算，最后求得近似极小值点(1.3175,-1.5675).满足误差ε=0.01.
2.1.5 多项式拟合
在反问题计算中，多项式拟合是常用的方法，其基本原理是：给定测量数据(xi,yi), i=1, 2, …,m, 求一个多项式
. y=a0+a1x+a2x 2+…+anx n
将数据代入得
上式可写为矩阵表达式：
Y=XA
这里：
两边乘以X 的转置X T 有
X T Y=X T XA
故有
编程：
clear all
t=[0,1,2,3] )
'()'(1Y X X X A -=n m
n m m m n n n n x a x a x a a y x a x a x a a y x a x a x a a y ++++=++++=++++=................................22102222210212121101⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=m n m m n n m a a a A x x x x x x X y y y Y .........1.........1...1 (21221121)
L=[0,9,35,90]
Y=L ’(MATLAB 中表示Y 的转置矩阵) for i=1:n
x(i,1)=1 for j=1:3
x(i,j+1)=x(i,j)*t(i) end end
A=inv(X ’*X)*X ’*Y
当n=1时，y=a+bx 为线性函数，可以由上式求出具体表达式：
式中E(X)为X 的平均值，E(Y)为Y 的平均值。

D(X)为X 的方差。

上式与最小二乘法得到的结果相同。

例. 已知数据(0,0), (1,1), (2,4),(3,8)， 求一元回归函数y=a+bx? 解. 我们求得m=4 E(X)=(0+1+2+3)/4=7/4;
D(X)=((-7/4)^2+(-3/4)^2+(1/4)^2+(5/4)^2)/4=1.3125 E(Y)=(0+1+4+8)/4=13/4
E(XY)=(0*0+1*1+2*4+3*8)/4=33/4 计算得： b=1.9524; a=-0.1667; 则一元回归为：
Y=-0.1667+1.9524x
2.1.6 数值积分
在数学模型竟赛中，能求出分析解的积分太少，大多只能用数值方法离散计算。

设h 为步长，
a=x0<x1<… <xn=b xk=a+kh zk=xk+0.5h
)
()()
()
()(X bE Y E a X mD Y E X mE y x
b k k
-=-=
∑
积分用复化中点公式计算
例.计算f(x)=x 2
在区间[-1,1]上的积分。

解.取步长h=1,得区间[-1,0],[0,1], 中点z0=-0.5, z1=0.5; 积分为
I=h[f(z0)+f(z1)]=(-0.5)^2+(0.5)^2=0.5
实际精确值为0.6667, 相差0.1667; 若取更小的步长，误差将变小。

对于二维积分，其积分区域为D ，将区域D 划分为若干四边形，第k 个四边形中点坐标为Xk, 面积为Sk, 则用复化中点计算
例. 计算f(x,y)=x 2
+y 2
在区域[0,2;0,2]]上的积分。

解.取步长h=1,得子区域[0,1;0,1],[0,1;1,2], [1,2;0,1],[1,2;1,2]; 中点为
P1=(0.5,0.5); P2=(0.5,1.5);P3=(1.5,0.5); P4=(1.5,1.5); 子区域面积Sk=1.0; 积分为
I=f(P1)S1+f(P2)S2+f(P3)S3+f(P4)S4 = f(P1) +f(P2)+f(P3)+f(P4)
=0.5^2+0.5^2+0.5^2+1.5^2+1.5^2+0.5^2+1.5^2+1.5^2 =10
实际精确值为10.667, 相差0.667, 相对误差为0.667/10.667=0.0625. 在可接受范围之内。

当然，取小的步长，能缩小误差值。

2.1.7 微分方程数值解
设初值问题 . y ' =f(x,y); y(x0)=y0 取步长为h, 点x1=x0+h, x2=x1+h, x3=x2+h, … 导数 y ’(x1)=(y1-y0)/h, 则有 . y1=y0+hf(x0,y0)
)}
(...)2()1()0({)()...()(11
21
10
--++++=+++=⎰⎰⎰⎰n xn
xn x x x x b a
z f z f z f z f h dx
x f dx x f ⎰⎰∑=D
k
k
S X
f dxdy y x f )(),(
. y2=y1+hf(x1,y1)
………
例. 给定初值问题, y’=1+xy, y(0)=1; h=0.1; 求y(0.2)? 解. y1=y0+hf(x0,y0)=1+0.1f(0,1)=1+0.1(1+0*1)=1.1
y2=y1+hf(x1,y1)=1+0.1f(0.1,1.1)=1+0.1(1+0.1*1.1)=1.111 则近似有y(0.2)=1.111
第三节反问题数学模型求解
例1 物体下落距离L与时间T，正问题是：已知物体的高度，测量下落时间，即t=t(x). 反问题是：已知物体下落时间，求物体的下落距离。

解设距离L与时间有多项式关系：
L=a+bT+cT2+dT3
选取取不同距离Li，测量下落时间Ti，i=1,2,…,n; 构造函数Q(a,b,c,d):
Q(a,b,c,d)=Σ[a+bTi+cTi2+dTi3-Li]2
取A=(a,b,c,d)T, 用最小二乘法，求出上式的极小值
A=（X T X）X T Y
这里Y=(L1,L2,…,Ln)T，X=(X1,X2,X3,X4);X1=(1,1,…,1) ', X2=(T1,T2,…,Tn) ',
X3=((T12,T22,…,Tn2) ', X4=((T13,T23,…,Tn3) '.
例如测量得一组数据为
表3.1 物体下落时间T与距离L
T L T L T L
0.1 0.04905 1.6 12.175 3.1 45.563
0.2 0.19571 1.7 13.739 3.2 48.545
0.3 0.43711 1.8 15.398 3.3 51.621
0.4 0.77304 1.9 17.152 3.4 54.792
0.5 1.2035 2 19 3.5 58.057
0.6 1.7284 2.1 20.942 3.6 61.417
0.7 2.3478 2.2 22.979 3.7 64.871
0.8 3.0618 2.3 25.11 3.8 68.42
0.9 3.8702 2.4 27.336 3.9 72.063
1 4.7731 2.5 29.656 4 75.802
1.1 5.7705
2.6 32.071 4.1 79.632
1.2 6.8624
2.7 34.581 4.2 8
3.563
1.3 8.0488
2.8 37.184 4.3 87.576
1.4 9.3298
2.9 39.883 4.4 91.712
1.5 10.705 3 4
2.676 4.5 95.874
由上面数据计算得：a=-0.00126; b=0.04469;c=4.7295;d=-0.0007893;
L=-0.00126+0.04469T+4.7295T2-0.0007893T3
上式近似于公式L=0.5gT2.
例2 年龄与身高。

正问题是，根据年龄T，每周岁测身高H，得到身高H与年龄T的关系H=H(T). 反问题是：已知身高H，求年龄T，即求关系T=H(T).
解设年龄T与身高H的关系为
T=a+bH
根据下面数据
表3.2 年龄T与身高H
T H T H
1 39 9 111
2 48 10 120
3 57 11 129
4 66 12 138
5 75 13 147
6 84 14 156
7 93 15 165
8 102 16 174
求得a=30, b=9; 即有
H=30+9T
例3 速度V与轨道形状z=f(x)，其摩擦系数为μ，z为高度，初始速度为V0，给定末速度为Ve=V(y=H),.求轨道形状z=f(x)。

解设y0=f(x0),y1=f(x1),在区间[x0,x1],物体速度近似取v0, 由能量守恒定理，设M为物体质量，V1=V(x=x1)为x=x1 处的速度，有
初始(动+势)能=结束(动+势)能+磨擦力作功Wf
其数学表达式为：
0.5MV02+gMy0=0.5MV12 +gMy1+Ef
Wf=μ(Mg cos (a1)+Rn)ΔS1 =μ[Mg(x1-x0)+Rn ΔS1]
式中a1为线段y0—y1与x 轴的夹角，Wf 为物体在x0--x1段移动时磨擦力作的功，ΔS1为线段的长度。

Rn 为离心力，定义为：
R 为曲率半径。

取y ’=(y1-y0)/h, y ”=(y2-2y1+y0)/h 2
. 将区间划分为[x0,x1],
[x1,x2],…,[xn-h,xn], xi=x0+ih, h 为步长，yi=f(xi)为对应函数值，当已知第i-1点的速度V i-1，则第i 点的速度Vi 为
0.5MVi 2
=0.5MV i-12
-gM(yi-y i-1)-Wf
Wf=μ(Mg cos (ai)+Rn)ΔSi=μ[Mg(xi-x i-1)+Rn ΔSi]
式中ai 为线段y i-1—yi 与x 轴的夹角，Ef 为物体在x i-1--xi 段移动时磨擦力作的功，ΔSi 为线段的长度。

将区间划分为[x0,x1],i=1,2,…,n. 当给定函数y=f(x)在点xi 时值yi, 用上面的公式能计算出末速度V(x=H). 设Y={y1,y2,…,yn}为y 在世点值组成的向量，与末速度V 有关系V=V(Y),给定Y=Yk, 得V=Vk, 则计算出数据(Yk,Vk), k=1,2,…,m, 令
Y=a+bV+cV 2
+dV 3
只要求出常数a,b,c,d, 对速度Ve, 我们就可以求出对应Y=Ye. 然后用插值方法求出曲线y=f(x).
考虑一个具体的例子，设y=f(x)定义在区间[0,1],且f(x=1)=0.5;末速度Ve=4, 初始速度V0=10, 摩擦系数为μ=0.01,单位取为焦尔，米，牛顿。

将[0,1]划分为两个区间[0,.5],[.5,1], 由于y0=0,y2=0.5, 则只有一个未知量y1. 改变y1 值，求得V 值。

2
/322
0)1(1y y K R R V M
Rn '+''=
==
表3.3 速度与坐标y1的关系
y1 v1(x=0.5) v2(x=1.0) 曲率半径R
0.025 4.9208 3.8725 0.68594
0.05 4.8679 3.9346 0.65706
0.075 4.8133 3.9947 0.61928
0.1 4.7565 4.0522 0.57137
0.125 4.6964 4.1065 0.512
0.15 4.6309 4.1561 0.43986
0.175 4.5549 4.1975 0.35365
0.2 4.4533 4.2208 0.2522
0.225 4.2481 4.1809 0.13454 设速度v=v2与坐标y1的关系
. y1=a+bv+cv2
我们求得a=7.6919; b=-4.2762; c=0.59305; 即有
. y1=7.6919-4.2762v+0.59305v2
取v=v2=4, 得y1=0.0759.
第四节具有反问题的数模竟赛题--单板滑雪
例4.1美赛2011问题A：单板滑雪
请设计一个单板滑雪场（现为“半管”或“U型池”）的形状，以便能使熟练的单板滑雪选手最大限度地产生垂直腾空。

“垂直腾空“是超出“半管”边缘以上的最大的垂直距离。

制定形状时要优化其他可能的要求，如：在空中产生最大的身体扭转。

在制定一个“实用”的场地时可能需要权衡哪些因素？
解单板滑雪场“U型池”见图4.1, 可以使单板爱好者们，从一面墙到另一面之间移动，跳跃并做花样动作。

墙面形状对运动速度有影响，运动轨迹是三维的，为了简化，只考虑二维截面。

由于对称性，只考虑一面墙，由于只有墙的形状可以改变，平底部分不能改变，因而墙面形状部分。

二维截面的墙面见图4.2, 图中L为墙面宽度，h’为高度，h为运动员离开墙面后上升高度。

G为运动员身体重量，Fu为离心力。

设墙面曲线为y=f(x), 运动员进入初始速度为V0, 末速度为Ve 。

当f(x)已知时，根据力学原理，能计算出末速度,这是正问题。

若给定末速度，求曲线f(x), 这是反问题。

可知，反问题比正问题难得多。

当曲线f(x)未知时，它要满足条件f(0)=0, f(L)=h ’, f ’(0)=0, 1/f ’(L)=0. 即运动员沿水平方向进入弧曲线，沿垂直方向离开弧曲线。

A1无磨擦
当不考虑磨擦损失时，由能量守恒原理，有
2
012mv mgH =
式中m 为运动员质量，H=h+h ’, 有H 的计算式：
2
02v H g =
故垂直腾空高度h=H-h ’. 它说明当不考虑磨擦时，垂直腾空高度与曲线形状无关。

A2有磨擦
当计及磨擦时，设运动员从点x=x0运动到点x=x1时，速度为v1，由能量守恒原理，有
Fnds
gh v v m μ=--)(5.0212
式中μ为磨擦系数，Fn 为作用在曲线法向上分力，ds 为曲线段长度,
θcos /)(1dx dx x f ds ='+=
且有
Fn=mgcos θ+mv 2
/R Fnds=mgdx+mv 2/Rds
θ为切线角，R 为曲率半径： 2
/32)1(1y y K R '+'
'=
=
取速度v=v0, 由上面公式计算得速度v1. 用类似的方法，将区间[0,L]划分为x0=0,x1,x2,…,xn=L. 依次计算出速度v1,v2,…,Vn. 再改娈曲线值
y1,y2,…,yn=h ’,求出对应速度，根据本题要求，要计算出最大速度Vn 对应的f(x)，则从所有计算结果中，选出最大的Vn 对应值y1,y2,…,yn 。

再用曲线拟合的方法计算f(x).
在实际计算中，μ=0.1,初始f(x)为圆弧
(
)(f x kL =
计算结果见表4.1, 由表可知，经过优化后，垂直腾空高度h 增加近20%.
表4.1: The comparison final velocity of initial shape and optimal
A3空气阻力
运动员速度高速运动时，将产生空气阻力：
22/1v kA R A '=ρ
K 为空气阻力系数，A 为迎风面积，ρ气体密度，v ’运动速度。

则有能量关系
A G F W W W E ++=∆
G F W W E ,,∆为损失动能、磨擦力作功、势能、阻力作功。

然后根据上面方法求出末速度Vn,和垂直腾空高度h 。

A4扭转
要求运动员腾空后，在空中产生最大的身体扭转。

产生的方法是：运动员用脚作用于U 池。

见图 4.3。

设E 为运动员离开U 型池作的功，E1为旋转能量，ω为旋转角速度， E2=E-E1为向上增加能量，E4为运动到池边，速度为Vn 的能量，E5=E4+E2为离开U 型池向上运动的能量
则有 则
在空中停留时间T 为
式中A 为运动员总转动角度。

将h 代入上式得
求得A 2
最大值时的ω值
再计算出V1，V2，V ，得角度为
Vn E m
V E mg h -=
=52
25
12
222
2
2
211)2(23661V V V V Vn m
E V m E L m E -=+===ωgh
A
T gh T 22==
=ω
a b 22=
ωE m
V b L a a b A n 2
3
/22
1422+
==-=ωωV
V 1)cos(=
α
计算结果见表4.2，由表可知，身体扭转角增加15%。

表4.2身体扭转
)
) )
Table 7: The improvement of setting slant angle
第五节具有反问题的数模竟赛题--储油罐
例5.1；2010年国赛A题储油罐的变位识别与罐容表标定
一内容
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，对于实际储油罐(见图5.1)，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b）之间的一般关系。

请利用罐体变位后在进/出油过程中的实际检测数据，根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm 的罐容表标定值。

进一步利用实际检测数据来分析检验你们模型的正确性与方法的可靠性。

二 逆问题数学模型
本文只考虑储油罐的纵倾和横倾变位，认为储油罐整体是钢性的。

设储油罐内液体体积V 与油浮显示油位高度h 之间的函数关系为
),,(βαh V V =
式中α为纵倾角，β为横倾角。

再设i V 为高度i h h =时的液体体积测量值,
1,2...i N =，定义理论计算值(,,)i V h αβ与测量值i V 的误差函数：
2]),,([),(i i V h V G -=∑βαβα (2.1)
由于环境条件，测量值对应的变位倾角00(,)αβ是未知的，但可以从上面公式求极值得到:
}]),,([{),(200i i V h V mim G -=∑βαβα (2.2)
当V 为倾角的线性函数时：
βαβα)()()(),,(210h a h a h a h V ++= (2.3)
上式中，()k a h ,0,1,2k =为已知函数，则有G 的表达式：
2210])()()([),(i i i i V h a h a h a G -++=∑βαβα (2.4)
利用最小二乘法，得极值解为：
2
00/)(/)(b
ad bg ah bh dg -=∆∆
-=∆-=βα (2.5)
式中：
)]
([)()]([)()
()()
()(0201212
221i i i i i i i i i i h a V h a h h a V h a g h a h a b h a d h a a -=-====∑∑∑∑∑ (2.6)
对于一般情形，当倾角较小时，取V 的线性项有较高的精确度，也能作为迭代初值。

由上面讨论可知，当已知罐内液体体积的测量值和理论计算值，变位倾角可以由极值计算。

下一节将讨论液体体积理论计算方法。

三．罐内液体体积理论计算
设Oxyz 为储油罐贴体坐标系，Oxy 平面为储油罐截面积最大的截面所在的平面，z 轴垂直于Oxy 平面。

再设液面所在的平面P 的法向量为n=(A,B,C)，它由倾角决定，对应平面方程为
)]()([1
0003y y B x x A C
z z -+--
= (3.1) 式中点(x 0,y 0,z 0)为油浮子在液面上的点。

本文假设油浮子始终在液面上。

平面Oxy 将储油罐截为上表面11(,)z f x y =, 和下表面22(,)z f x y =。

设Dxy 为平面z=0在储油罐内部区域，则罐内液体体积V 为：
⎰⎰=
Dxy
dxdy z z z F h V ),,(),,(3
2
1
βα (3.2)
根据平面P 与油罐上下表面之间的关系，式中F 由下面公式表示出：
⎪⎩
⎪
⎨⎧<>-≥≥-=2313212312
33210.....),,(z z z z z z z z z z z z z z F (3.3)
由于罐内液体体积V 为二重积分，当油罐上下表面用简单函数表示时，我们能求出积分表达式。

当罐面比较复杂时，可以用数值方法计算。

在下面实例中，可将平面Dxy 以步长dxdy 分成N 个小区域，将每一个区域对应的油柱近似为直方体，然后对所有N 个直方体进行离散求和。

四 卧式圆筒球形封头型油罐
4.1 上下表面和油面
卧式圆筒球形封头型油罐的空间直角坐标系见图1, 图中取Oxyz 为贴体坐标系，中截面在Oxy 平面上，取为上节中二重积分区域Dxy, L 为储油罐主体（圆柱体）的长度，r 为储油罐主体的截面园半径，q 为两端球冠体的高。

图1 油罐的空间直角坐标系的建立
设油浮子始终作为油面上的点，坐标为 (,0,)2
L
l h r -+-，这里h 为油浮显示油位高度，储油罐纵倾角为α，横倾角为β。

4.2 计算结果
根据前面定义的误差函数(,)G αβ, 利用2010全国大学生数学建模竞赛A 题(/)所给出的实际测量的前302个数据，采用变步长遍历搜索得误差函数(,)G αβ最小值对应的油罐纵倾角为α=2.1度，横倾角为β=4.3度. 利用部分测量数据计算参数，另一部分测量数据检验计算的参数是科学研究中常用的方法[8]。

4.3检验变位参数,αβ正确性
运用求得的α与β的值得到体积V 与h 的对应关系，对A 题所给出的实际测量剩余的后299个数据进行检验。

则模型所得数据与实际数据的相对误差，如下图10所示。

00.51
1.5
2 2.53
-0.05
-0.04-0.03-0.02-0.0100.010.020.03
0.04出油量的相对误差
显示油高（m ）
图2 模型计算数据与实验数据相对误差
由图中相对误差图可知，模型的计算数据与实验数据误差不超过5%，或计算值与实测值吻合精度高于95%. 实际计算得标准偏差为0.0348, 平均绝对误差小于0.01,可以认为所建立的模型是正确的。

4.4 罐容表标定值
油罐变位后必需重新标定，根据本文计算公式，计算的标定值见表1, 油位高度间隔取为10cm 。

表1 变位后油位高度间隔为10cm 的罐容表标定值
注：表格中油面高度的单位为米（m），油的容量的单位为升（L）五结论
根据油面与油罐表面之间的相对关系，选用油罐贴体坐标系，求出油量V与油面高度及变位倾角之间的计算关系式，然后利用测量数据，构造逆问题求解倾角的误差函数。

求出了极小值对应的倾角值。

最后对罐容表进行了重新标定。

文中油量V为二重积分值，有多种数值方法能给出满足实际要求精度的计算值，且与用解析积分表达式计算值接近。

当油罐表面函数比较简单时，可以求出积分表达式，或部分区域上的积分表达式。

但当积分表达式较繁杂时，编写计算程序难度提高。

因而，用数值方法计算积分是一个可行的方法。

文中方法不仅适用于卧式圆筒球形封头型油罐，也能用于其它类型油罐逆问题计算。

第六节 2012B题太阳能小屋的设计
1、问题的重述
太阳能小屋，就是在小屋的顶面及外墙铺设光伏电池，通过光伏电池吸收太
阳辐射能量转化为直流电，在经过逆变器转化为 220 的通用交流电源供居民生活使用，多余的存入电网。

光伏电池对太阳的吸收与电池的种类，辐射强度，阳光直射角度等诸多条件
有关，而不同种类的光伏电池的价格也相差甚远。

根据附件中的数据，在考虑年
发电总量又考虑铺设成本的情况下对小屋电池排布进行设计，计算出小屋光伏电
池 35 年寿命期内的发电总量、经济效益（当前民用电价按 0.5 元/kWh 计算）及投资的回收年限，并选出合适的逆配器与串并联接线方式。

问题 1：根据山西省大同市的气象数据，仅考虑贴附安装方式，选定光伏电
池组件，对小屋的部分外表面进行铺设，并根据电池组件分组数量和容量，选配
相应的逆变器的容量和数量。

问题 2：电池板的朝向与倾角均会影响到光伏电池的工作效率，请选择架空
方式安装光伏电池，重新考虑问题 1。

问题 3：根据给出的小屋建筑要求，请为大同市重新设计一个小屋，要求画
出小屋的外形图，并对所设计小屋的外表面优化铺设光伏电池，给出铺设及分组
连接方式，选配逆变器，计算相应结果。

问题一的分析
我们在设计电池排布的时候既要考虑太阳辐射的能量，又要考虑铺设成本。

然而题中已知每千瓦时的电价，我们可以将设备 35 年寿命内转换的电量转化为收益，然后减去初始成本得到总收益，作为目标函数简化求解。

在对每一块墙面上电池板的排布，由于电池板的尺寸，大小和形状皆不经相同，这样一个排布上的优化问题就涉及到排样算法，十分复杂而且只能不适合于本题的目标函数。

于是我们考虑将矩形拼排的问题转化为长度问题。

经过观察发现大部分电池板宽度都在一米左右浮动，于是我们将墙面进行切割，结合实际可用空间形状切割为宽度一米左右的若干个条状得到总长度，然后要求所有电池板的长度和的数值小于该总长度即可。

为保证光伏组件正常工作，考虑到光伏电源串并联的要求：只允许相同型号的光伏组件进行串联。

多个光伏组件串联后可以再进行并联，并联的光伏组件端电压相差不应超过 10%。

最后连接，形成组件阵列。

最后考虑到逆变器的额定工作电压（V）范围和功率容量（W）等参数进行分组设计。

问题一的求解：
首先我们要从房屋的六个表面中选出几个铺设电池板的墙面，而且不考虑四
周其他建筑的遮光效应，如果一个面被选定则该面整面都可以用来铺设电池板。

我们将六个面编号。

图6-1太阳能小屋
六个面单位时间接收到的辐射等于每个面在水平面、东、南、西、北五个面
投影面积乘以该方向的辐射强度。

首先不考虑门窗，根据题目图中数据，六个面在五个方向的投影面积为：
表 1 墙面投影面积表（单位：平方米）
根据题中给出的太阳辐射强度，计算全年每个面接收到的能量。

在铺设同等面积电池板时，求得全年能接受到的太阳辐射总能量墙面②最多，而墙面⑤最少，为了保障铺设电池板的经济效益，我们优先考虑按照②、①、⑥、④、③、⑤顺序安排电池板。

然后我们从墙面②开始，由于实际电池块大部分小于1.114，我们先对墙面
进行划分
将区域分为 6 个条状长方形，得到几个长方形长的总和，使用的所有的电池矩形的长的数量和应该小于这六个长方形长的数量和。

最后结果为：。