姜启源 谢金星《数学建模案例集》.pdf

【案例1】原子弹爆炸的能量估计1问题的提出1945年7月16日，美国科学家在新墨西哥州阿拉莫戈多沙漠(L0 s Alam)进行了“三位一体实验”(Trinity Tes),试爆了全球第一颗原子弹（图1）.这一事件令世界为之震惊，井从某种程度上改变了第二次世界大战以及战后世界的历史但在当时，有关原子弹爆神的任何货料都是保密的，一般人无法得到任何有关的数据或影像货料，因此人们无法比较准确地了解这次燥炸的威力究竞有多大两年以后，美国政府首次公开了这次爆炸的录影带，但没有发布任何有关的数据.英国物理学家G.L.Tayo(1886一1975)通过研究这次爆炸的录影带，建立数学模型对这次爆炸所释放的能量进行了估计，得到的估计值为19,2千吨（干吨即相当于1千吨TNT的核子能量)，后来正式公布的信息显示，这次爆炸所释放的实际能量为2】千吨，可见两者是相当接近的，雷1原子弹浮炸示意图部么，Taylor是如何对原子弹爆炸的能量进行估计的呢？Taylor知道，爆均产生的冲击被以爆炸点为中心呈球面向四周传插，爆炸的能量越大，在一定时刻冲击波传播得越运，而冲击被又可以通过爆炸形成的“警酷云”反陕出来.Tyo研究这次爆炸的录影带，测量出了从爆炸开始，不同时刻爆炸所产生的“曹菇1云”的半径大小.表1是他测量出的时刻t所对应的“蘑菇云”的半径r(t).现在的任务就是利用表1和其他知识，估计这次爆炸所释放的能量表1时刻(ms)所对应的“蘑菇云”的半径r(m】r()tr()r(t)r(t)r(t)0.1011.10.8034.21.5044.43.5361.115.0106.50.2419.90.9436.31.6546.03.8062.925.0130.00.3825.41.0838.91.7946.94.0764.334.0145.00.5228.81.2241.01.9348.74.3465.653.0175.00.6631.91.3642.83.2659.04.6167.362.0185.0建立数学模型表1只给出了爆炸后不同时刻“蘑菇云”的半径，由此无法直接计算爆炸能量，必须建立模型将不同时刻“蘑菇云”的半径与爆炸能量联系起来Taylor建立计算爆炸能量的数学模型所采用的是量纲分析法，记爆炸能量为E,将“蘑菇云”近似看成一个球形形状.除时刻：和能量E外，与“蘑菇云”的半径r有关的物理量还可能有“蘑菇云”周围的空气密度（记为）和大气压强(记为P),将要寻求的关系r=(t,E,p,P)(1)记作更一般的形式f(r,t,E,p,P)=0,(2)其中有5个物理量，接下来的任务是用量纲分析法确定这个函数关系2取3个基本量纲：长度L,质量M和时间T,(2)中各个物理量的量纲分别是r=L,t=T,E=LMT2,p=L-3M,P=LMT-2,(3)由此得到量纲矩阵r102-3-1A35001(4)L01-20-2因为A的秩是3，所以齐次方程2Ay=0,y=(y,y,y3,y,3)*(5)有5-3=2个基本解令y1=J,y3=0,得到一个基本解y=(1,-2/5,-1/5,1/5,0)；令y,=0,y=1,得到另一个基本解y=(0,6/5,-2/5,-3/5,1)T.根据量纲分析的Bucking-ham Pi定理，由这2个基本解可以得到2个无量纲量m,=gp=r(品)15(6)、1/5T2=t5E-2p35P=Eo(7)且存在某个函数F使得F(T1,T2）=0(8)与(2)等价.为了得到形如(1)的关系，取(8)的特殊形式1=中(2)（其中中是某个函数)，由(6)、(7)即(9)于是(10)函数中的具体形式需要采用其他方式确定.(10)式就是量纲分析法确定的、如(1)式所要求的函数关系.3.数值计算为了利用表1中：和r的数据由(10)确定由于原子弹爆炸的能量E,必须先估计无量纲量也（，）的大小Taylor认为，对于原子弹爆炸来说，所经历的时间非常短，而所释放的能量是非常大的.仔细分析(7)可知，m=(。】0.他又根据一些小型爆炸试验的数据，最终建议取(0)1，这就是说，(10)可以近似为I/5P)(11)这一近似关系式表明，半径与大气压强P无关，而当E,P一定时r与2”成正比.我们可以先用表1的数据检验一下这个关系，设36.95+6.96.85o1-6.86.75斜6.76.656.6-4-3.5-2.51.5logof图3(14)式与实际数据(+)的拟合(为y=6.9038)对它作最小二乘拟合相当于取(15)右端的平均值，得到E=8.282510,爆炸能量是19.7957（千吨）.从这个例子可以看出，Taylor采用的量纲分析方法在这里取得了巨大成功.附录MATLAB程序clear,clc;data=0.1011.10.8034.21.5044.43.5361.115.0106.50.2419.90.9436.31.6546.03.8062.925.0130.00.3825.41.0838.91.7946.94.0764.334.0145.00.5228.81.2241.01.9348.74.3465.653.0175.00.6631.91.3642.83.2659.04.6167.362.0185.0;t=r=;for i=1:5t=t:data(:,2*i-1);%得到t的列向量r=r;data(:,2*i);%得到r的列向量endrl=log(r);t1=log(t);s polyfit(tl,r1,1)%(12)取对数后的线性最小二乘拟合5b=s(1),%（12)的系数ba=s(2);t=0:0.1:70;rr exp(a)tt.b;%用拟合结果由(12)进行计算plot(tt,rr,t,r,+)pausex=1og10(t*1e-3);%t的单位由ms化为$y=5/2*og10(r）-x;plot(x,y,+)xlabel(log10(t);ylabel(5/2 log10(r)-log10(t);c=mean(y)）%y的平均值hold on;plot(x,e,.-);hold off;rou0=1.25;E=0u0*10(2*c)kiloton=E/4.184e12E1=rou0*mean(r.5./L.2)*1e6%用(15)作拟合kiloton1 E1/4.184el2参考文献1Illner R,Bohun C S,Mc Collum S,Van Roode T.Mathematical Modelling:ACase Studies Approach.Student Mathematical Library,AMS(American Mathe-matical Society)2005,27.2姜启源，谢金星，叶俊，数学模型.第三版.北京：高等教育出版社，2003.16