河南城建学院

0.0150.010.0050-0.005-0.01-0.01501234567

(3)大气压强p随高度x变化的理论公式为

，为验证这一公式，

测得某地大气压强随高度变化的一组数据如表所示。试用插值法和拟合法进行计算并绘图，看那种方法较为合理，且总误差最小。

高度/m 压强/Pa 0 0.9689 300 0.9322 600 0.8969 1000 0.8519 1500 0.7989 2000 0.7491 插值法： >> x=[0 300 600 1000 1500 2000]; p=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; >> xi=linspace(0,2000);p0=1.0332*exp(-(xi+500)/7756); >> p1=interp1(x,p,xi,'spline'); >> plot(xi,p0,'*',xi,p1) >> e1=p1-p0; >> e=sum(e1.^2) e =

1.8652e-005 拟合法：

>> x=[0 300 600 1000 1500 2000]; p=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; >> P=log10(p)

第 9 页

河南城建学院

P =

-0.0137 -0.0305 -0.0473 -0.0696 -0.0975 -0.1255 >> p1=polyfit(x,P,1) p1 =

-0.0001 -0.0137

>> b=p1(1)/0.4343,a=10.^p1(2) b =

-1.2863e-004 a = 0.9689

>> xi=linspace(0,2000);p0=1.0332*exp(-(xi+500)/7756); >> p2=polyval(p1,xi);P2=10.^p2; >> e2=P2-p0;e=sum(e2.^2) e =

1.8116e-005

第 10 页

河南城建学院

四 Matlab数值运算

1、数值微积分

练习：瑞士地图如图所示，为了算出其国土面积，首先对地图作如下测量：以由西向东方向为X轴，由南到北方向为Y轴，选择方便的原点，并将从最西边界点到最东边界点在X轴上的区间适当划分为若干段，在每个分点的Y方向测出南边界点和北边界点的Y坐标Y1和Y2，根据地图比例尺知道18mm相当于40km，试由测量数据计算瑞士国土近似面积，与其精确值41228km比较。 X Y1 Y2 X Y1 Y2 >>

x=[7,10.5,13,17.5,34,40.5,44.5,48,56,61,68.5,76.5,80.5,91,96,101,104,106.5,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,116,122,83,81,82,86,85,68];

>> X=x./18*40;Y1=y1./18*40;Y2=y2./18*40; >> t1=trapz(X,Y1),t2=trapz(X,Y2), t=t2-t1 t1 =

3.3819e+004 t2 =

7.6328e+004 t =

第 11 页

96 43 101 104 106.5 111.5 118 123.5 136.5 142 146 150 157 158 37 33 28 32 121 65 55 54 83 52 81 50 82 66 86 66 85 68 68 7 44 44 10.5 13 45 59 47 70 17.5 34 50 72 50 93 40.5 44.5 48 38 30 30 56 34 61 36 68.5 76.5 80.5 91 34 41 45 46 2

100 110 110 110 117 118 116 118 118 121 124 121 121 116 122 河南城建学院

4.2510e+004 >> expt=t-41228 expt = 1.2819e+003 2、习题

(4)利用梯形法和辛普森法求定积分的值，并对结果进行比较。如果积分区间改为-5~5结果有何不同？梯形积分中改变自变量x的维数，结果有何不同？

>> x=linspace(-3,3);y=exp(-x.^2/2); >> t=(1/2*pi)*trapz(x,y) t = 3.9267 >>

q=(1/2*pi)*quad('exp(-x.^2/2)',-3,3)

>> x=linspace(-5,5);y=exp(-x.^2/2); >> t=(1/2*pi)*trapz(x,y) t = 3.9374 >>

q=(1/2*pi)*quad('exp(-x.^2/2)',-5,5)

>> x=linspace(-3,3,150);y=exp(-x.^2/2); >> t=(1/2*pi)*trapz(x,y) t = 3.9268

(5)分别用矩形法、梯形法、辛普森法和牛顿-科茨4种方法近似计算定积分取n=4，保留4位有效数字。

矩形法：

>> x=linspace(0,1);y=x./(x.^2+4);

>> t=cumsum(y)*1/99;T=t(100) T =

，

q = 3.9374 q = 3.9268

第 12 页