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

数学模型中的反问题向下运动向上运动风筝数学模型竟赛中有很多涉及反问题。如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。如用一维传热公式：反问题是：已知温度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|L1, 则迭代收敛。证：由于x*=g(x*)，则xk+1-x*=g(xk)-g(x*)=g() (xk-x*)有| xk+1-x*|L| xk-x*|L2| xk-1-x*|Lk+1| x0-x*|因为L0, 故xkx*. 证毕。例 求f(x)=x-x*x的零点。解：这里g(x)=x*x, g(x)=2x, 则当|x|0.5时，|g(x)|1, 即|x|0.5时，迭代公式. xk+1=xk2收敛。取x0=0.1, 计算得X1=x02=0.12=10-2X2=x12=(10-2)2=10-4.最后求得xkx*=0. 实际上，我们知道x=0为x=x*x的解，但它还有一解x=1; 由于|2x|=|2*1|=2, 则用上面迭代公式x=g(x)=x*x求不出解x=1. 它需要构造另一种迭代公式.xk+1=g(xk)=xk容易验证当x=1时，|g|1. 取x0=2, 计算得X1=x00.5=20.51.414X2=x10.5=(20.5) 0.5=20.251.1189X3=x30.5=(20.25) 0.5=20.1251.090.最后求得xkx*=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, xix*. 例. 求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, 则零点x*在区间x2,b1=0.875,1.25中，故取a2=0.875, b2=1.25. 如此计算下去， 当bi-ai0, 无根，停止计算。否则转下一步；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; 否则a1a0, b1b0, 转第三步。二分法能用图形来说明，其示意图见图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)F(d0-h), 取d1=d0-h; 否则取d1=d0+h;5)令d1d0, h2h, 转第2步。6)取a=d0-h, b=d0+h, 用二分法求极值点。二分法求极值点的原理与求根原理类似。由下面定理给出：定理：设F(x)在a,b上连续，且0F(x),若c为a,b中的点，且 F(c)minF(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)=minF(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)=minF(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. 可以证明，xix*为极小值点。由上面讨论可知，求极小值点分为两步，先求极点所在的区间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=minF(X0),F(X-h),F(X+h), 和对应的坐标点X1, 若X1=X0，取hh/2, 步长减半，否则取h2h, 步长加倍，再将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)-12+y(x+y)-y-22取初值点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+a2x2+anxn将数据代入得上式可写为矩阵表达式： Y=XA这里：两边乘以X的转置XT有 XTY=XTXA故有编程：clear allt=0,1,2,3L=0,9,35,90Y=L(MATLAB中表示Y的转置矩阵)for i=1:nx(i,1)=1for j=1:3 x(i,j+1)=x(i,j)*t(i)endendA=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=4E(X)=(0+1+2+3)/4=7/4; D(X)=(-7/4)2+(-3/4)2+(1/4)2+(5/4)2)/4=1.3125E(Y)=(0+1+4+8)/4=13/4E(XY)=(0*0+1*1+2*4+3*8)/4=33/4计算得： b=1.9524; a=-0.1667; 则一元回归为：Y=-0.1667+1.9524x2.1.6 数值积分在数学模型竟赛中，能求出分析解的积分太少，大多只能用数值方法离散计算。设h为步长，a=x0x1 New-M-file3)在编辑窗中输入程序；4）点击Debug-Save and Run5）输入文件名w11, 保存6）屏幕上出现结果c=3;上机实习程序二；计算 S=1+2+.+10clear all s=0;for k=1:10; s=s+k;end;sun=s上机实习程序三计算s=2+4+6+.+20; 如s10, 输出“优秀”，否则输出“下一次优秀”clear alls=0;for k=1:10; s=s+2*k;end;if(s10);display(优秀),end if(s4);continue;end; f2=x1+4*x2+x3; if(f14);continue;end; f3=2*x2+x3; if(f36);continue;end; f4=x1+x2; if(f41);continue;end; z=3*x1-x2+4*x3; if(Mzz);Mz=z;y(1)=x1;y(2)=x2;y(3)=x3;end; end;%i3 end;%i2end;%i1y(4)=Mz程序中continue为进行下一个循环，而不计算后面的语句。程序二：clear allf=-3;-1;4;A=1 3 -2 1 4 1 0 2 1 1 1 0;b=4;4;6;1;x，fval=bintprog(f,A,b)上面bintprogo为0-1整数线性规划函数，求解的为目标函数的最小值, 所以要在f前面加个负号。clear allu(1)=1.3175;u(2)=-1.555;h=0.0125;k=0; for j=1:2;z=u;for i=-1:1;z(j)=u(j)+h*i; x=z(1);y=z(2); k=k+1;a(k,1)=x;a(k,2)=y; a(k,3)=(x-y*(x+y)-1)2+(y*(x+y)-y-2)2; end; end; aclear allv0=5;g=9.81;M=1;mu=0.01;y1=0.25;y0=0;y2=0.5;h=0.5; x0=0;x1=0.5;x2=1;for k=1:10;y1=0.025*k; yd1=(y1-y0)/h;yd2=(y2-2*y1+y0)/h/h; r=yd2/(1+yd12)1.5; ds=sqrt(x1-x0)2+(y1-y0)2)0.5;w=0.5*v02-g*y1-mu*(g*0.5+v02/r*ds);if(w0);break;end;v1=sqrt(2*w); yd1=(y2-y1)/h;yd2=(y2-2*y1+y0)/h/h; r=yd2/(1+yd12)1.5; ds=sqrt(x2-x1)2+(y2-y1)2)0.5;w=0.5*v02-g*(y2-y1)-mu*(g*0.5+v12/r*ds);if(w0);break;end;v2=sqrt(2*w); v(k,1)=y1;v(k,2)=v1;v(k,3)=v2;v(k,4)=r;end; vfor i=1:9;w=v(i,3);x(i,1)=1;x(i,2)=w;x(i,3)=w*w;end;b=x*x; b1=x*v(:,1);bv=inv(b);b2=bv*b136