数值积分与数值微分编程计算.docVIP

解：卫星轨道的示意图如右上图所示，分别是椭圆轨道的长半轴和短半轴，地球位于椭圆的一个焦点处，焦距为，地球半径为，近地点和远地点与地球表面的距离分别是和. 由图中可知，上述数据存在如下关系： 由椭圆性质 ，将的数据代入以上各式可得，. 椭圆的参数方程为： ， 根据计算参数方程弧长的公式，椭圆长度可表为如下积分： 由于该积分无法求得解析解，下面我们编写MATLAB程序对其进行数值求解。 s1=439;s2=2384;r=6371; a=(s1+s2)/2+r a = 7.7825e+003 c=a-s1-r; b=sqrt(a^2-c^2) b = 7.7215e+003 y=inline(sqrt(7782.5^2*sin(t).^2+7721.5^2*cos(t).^2)); %建立积分内联函数 t=0:pi/10:pi/2; y1=y(t); format long L1=4*trapz(t,y1) %梯形积分 L1 = 4.870744099902405e+004 L2=4*quad(y,0,pi/2,1e-6) %辛普森积分 L2 = 4.870744099903280e+004 求解结果显示：两种方法计算求得的积分结果相当接近，轨道长度约为：. 解：我们需要求出上图中不规则图形的面积，而根据积分的定义可知这实际上就可以归结为一个积分问题，我们采用梯形公式对其进行数值积分，MATLAB程序代码如下： x=[7.0 10.5 13.0 17.5 34 40.5 44.5 48 56 61 68.5 76.5 80.5 91 96 101 104 106 111.5 118 123.5 136.5 142 146 150 157 158]; y1=[44 45 47 50 50 38 30 30 34 36 34 41 45 46 43 37 33 28 32 65 55 54 52 50 66 66 68]; y2=[44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85 68]; y=y2-y1; format long s=trapz(x,y)*40^2/18^2 %梯形积分 s = 4.241481481481482e+004 结果表明：瑞士国土约为.(本题也可采用其它积分方法计算面积，如分段线性插值，辛普森积分等). 解：（1）z=2*quad(exp(-x.^2)./(1+x.^4),0,100000,1e-6) z = 43484833213566 (2) z=quad(sin(x)./(1-x.^2).^(1/2),0,0.99999,1e-6) z = 88948175020513 (3) z=quad(1./(x.^0.5.*(1+sin(x))),1e-6,1,1e-6) z = 1.58462649585356 （4）本题可将对原积分参数进行变换，令，，则原积分可化为二次积分：，然后用dblquad命令进行求解，亦可采用蒙特卡罗方法直接进行数值积分,下面分别用这两种方法进行求解. (I) f=inline((2+r*cos(a)+r*sin(a))*r); z=dblquad(f,0,2*pi,0,1) z = 6.28318531935223 (II) n=100000;u=0;m=0; x=unifrnd(-1,1,1,n); y=unifrnd(0,2,1,n); for i=1:n if x(i)^2+y(i)^2=2*y(i) u=u+1+x(i)+y(i); end end p=4*u/n p = 6.27878234230698 事实上，对于这一积分，我们可以求得其精确值为. 从上面计算的结果容易看到，方法（I）的精度很高；方法（II）的精度较差，但其优点在于不需要对原积分进行转换. 解：先计算蛋糕的体积 s=inline(pi.*(2-cosh(2.*h).*2/5).^2); v=quad(s,0,0.999999999) v = 5.4171 由于题中未给出蛋糕实际上是怎样配料的，所以剩下的工作无法进行. 当然我们也可以自己设定有关数据进行计算，这里就不再赘述了.