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

解：卫星轨道的示意图如右上图所示，,a b 分别是椭圆轨道的长半轴和短半轴，地球位于椭圆的一个焦点处，焦距为c ，地球半径为r ，近地点和远地点与地球表面的距离分别是1s 和2s . 由图中可知，上述数据存在如下关系：

12122,a s s r c a s r =++=--

由椭圆性质 b =，将12,,s s r 的数据代入以上各式可得7782.5a km =，7721.5b km =.

椭圆的参数方程为：

c o s ,s i n x a t y b t == ，

(02)t π≤<

根据计算参数方程弧长的公式，椭圆长度可表为如下积分：

/2

22221/20

4(sin cos )L a t b t dt π=+⎰

由于该积分无法求得解析解，下面我们编写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

求解结果显示：两种方法计算求得的积分结果相当接近，轨道长度约为：4

4.8710km ⨯.

解：我们需要求出上图中不规则图形的面积，而根据积分的定义可知这实际上就可以归结为一个积分问题，我们采用梯形公式对其进行数值积分，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

结果表明：瑞士国土约为42

4.241510km ⨯.(本题也可采用其它积分方法计算面积，如分段线性插值，辛普森积分等).

解：（1）z=2*quad('exp(-x.^2)./(1+x.^4)',0,100000,1e-6)

z =

1. 43484833213566

(2) z=quad('sin(x)./(1-x.^2).^(1/2)',0,0.99999,1e-6)

z =

0. 88948175020513

(3) z=quad('1./(x.^0.5.*(1+sin(x)))',1e-6,1,1e-6)

z =

1.58462649585356

（4）本题可将对原积分参数进行变换，令cos ,1sin x t y t ==+，(02)t π≤<，则原积分

可化为二次积分：21

00

(2cos sin )dt r t r t rdr π++⎰⎰，然后用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

事实上，对于这一积分，我们可以求得其精确值为2 . 从上面计算的结果容易看到，方法（I）的精度很高；方法（II）的精度较差，但其优点在于不需要对原积分进行转换.

解：先计算蛋糕的体积

s=inline('pi.*(2-cosh(2.*h).*2/5).^2');

v=quad(s,0,0.999999999)

v =

5.4171

由于题中未给出蛋糕实际上是怎样配料的，所以剩下的工作无法进行. 当然我们也可以自己设定有关数据进行计算，这里就不再赘述了.