梯形法数值积分

《MATLAB 程序设计实践》课程考核

一、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精

通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）

“矩形法、梯形法数值积分”

1.梯形法数值积分

A．算法说明：

梯形法数值积分采用的梯形公式是最简单的数值积分公式，函数梯形法数值积分表达式为：

f(x)在区间[a,b]上计算

baf(x)dxba[fa()fb( )]2由于用梯形公式来求积分十分粗糙，误差也比较大，后来改进后提出了复合梯形公式：

hban,其中，n为积分区间划分的个数；h为积分步长。

在MATLAB中编程实现的复合梯形公式的函数为：Combine Traprl. 功能：复合梯形公式求函数的数值积分。 调用格式：[I,step]=CombineTraprl(f,a,b,eps). 其中，f为函数名； a为积分下限； b为积分上限； eps为积分精度； I为积分值；

Step为积分划分的区间个数

B．流程图

开始 n=1; h=(b-a)/2; I1=0; I2= ( subs (sym ( if ) , findsym ( sym ( f ) ) , a ) + subs ( sym ( f ), findsym ( sym ( f ) ) , b ) ) / h; 是 abs ( –I1 ) >eps I2 否 否 n=n+1; h=(b-a) / n; I1=I2; I2=0; for i=0:n-1 是 x=a+h*i; x1=x+h; I2=I2+( h/2 ) * ( subs ( sym ( f ) , findsym ( sym ( f ) ) , x ) +subs (sym ( f ) , findsym ( sym ( f ) ), x1) ); I=I2; Step=n; 结束 C．复合梯形公式的MATLAB代码：

function[I,step]=CombineTraprl(f,a,b,eps) %复合梯形公式求函数f在区间[a,b]上的定积分 %函数名：f %积分下限：a %积分上限：b %积分精度：eps %积分值：I

%积分划分的子区间个数：step

if(nargin==3)

eps=1.0e-4; %默认精度为0.0001 end n=1; h=(b-a)/2; I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0;

for i=0:n-1 %第n次的复合梯形公式积分

x=a+h*i; %i=0 和n-1时，分别代表积分区间的左右端点 x1=x+h;

I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n;

D．应用举例:

复合梯形法求数值积分应用举例，利用复合梯形法计算定积分a. 流程图

开始 421dx 2x1f= 1/(x^2-1) a=2; b=4 调用函数CombineTraprl(f,a,b,eps) 结束

b. 原程序代码：

[q,s]=CombineTraprl('1/(x^2-1)',2,4) %精度为默认的10-4 结果

[q,s]=CombineTraprl('1/(x^2-1)',2,4,1.0e-6) %精度为10-6 结果

所以从复合梯形公式可以得出

421dx0.2939 2x12．矩形法数值积分的源程序

function[I,step]=CombineTraprl(f,a,b,eps) % 复合矩形公式求函数f在区间[a,b]上的定积分 %函数名：f %积分下限：a %积分上限：b %积分精度：eps %积分值：I

%积分划分的子区间个数：step if(nargin==3)

eps=1.0e-4; %默认精度为0.0001 end n=1; h=b-a; I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0;

for i=0:n-1 %第年n次的复合矩形公式积分

x=a+h*i; %i=0 和n-1时，分别代表积分区间的左右端点 x1=x+h;

I2=I2+h*subs(sym(f),findsym(sym(f)),x1); end end I=I2; step=n;

应用举列：复合矩形法求数值积分应用举例，利用复合矩形法计算定积分

21x2dx

a.流程图

开始 f= x^2 a=1; b=2 调用函数CombineTraprl2(f,a,b,eps) 结束

b.原程序代码：

[q,s]=CombineTraprl2('x^2',1,2) %精度为默认的10-4 结果：

二、科学计算和工程实际问题和举例

1.（题目）将100个学生5门功课的成绩存入矩阵P中，进行如下处理： （1）分别求每门课的最高分、最低分及相应学生序号。 （2）分别求每门课的平均分和标准方差。

（3）5门课总分的最高分、最低分及相应学生序号。

（4）将5门课总分按从大到小顺序存入zcj中，相应学生序号存入xsxh。 流程图：

开始 构造一个矩阵做学生序号 再构造一个随机矩阵做学生成绩 让每门课最高分放入a矩阵，最低分放人b矩阵 使用函数max求出每门课最大值，并用find函数找到，然后输入到a；同理用min求出每门课最小值。 使用mean函数求出每门课平均分c 使用var函数求出每门课的标准方差 使用sum函数求和，并找到每门课的最高分放入a2矩阵,每门课最低分放入b2。 最后用sort函数和descend把总分排序放人zcj文件和xsxh文件 结束

源程序代码：

clear

num=[1:1:100]';

p=[num,floor(50*rand(100,5)+45)]; a=[];b=[]; for i=1:5

m=max(p(:,i+1));

pp=find(p(:,i+1)==max(p(:,i+1))); aa=ones(size(pp))*m; a=[a;0,0;pp,aa]; end for i=1:5

m=min(p(:,i+1));

pp=find(p(:,i+1)==min(p(:,i+1))); bb=ones(size(pp))*m; b=[b;0,0;pp,bb];

end

c=mean(p(:,2:6)); d=var(p(:,2:6)); ps=sum((p(:,2:6))'); ps=ps';

a2=[find(ps==max(ps)),max(ps)]; b2=[find(ps==min(ps)),min(ps)]; a b c d a2 b2

[zcj,xsxh]=sort(ps,'descend')

结果：

（1）每门课最高分、最低分及相应序号。

（2）每门课的平均分和标准方差

（3）5门课总分总高分、最低分及相应学生序号

（4） zcj =

400 397 396 389 389 386 384 384 382 379 378 377 377 376 375 374 374 373 371 371 368 367 366 365 365 364 364 364 363 362 362

362 361 361 360 358 357 357 355 353 352 352 352 352 351 351 351 350 349 349 348 348 348 348 347 345 345 344 343 342 342 341 341 340 339 338 338 336 335 335 335 333 330 330 330

329 328 328 327 325 323 322 322 317 315 315 314 314 312 311 310 310 308 303 299 287 285 277 275 267

xsxh =

85 93 68 15 97 26 18 53 89 95 39 14 83 9 56

13 63 77 24 34 44 30 51 2 7 62 82 98 50 59 70 86 10 79 22 64 87 88 43 25 8 16 32 92 20 41 58 91 12 31 4 37 46 73 67 42 57 84 6

47 72 5 100 80 66 48 54 75 27 29 55 23 17 45 49 76 1 71 78 38 21 11 52 3 36 94 28 81 90 61 74 99 33 65 60 19 69 40 96 35

2．某气象观测站测得某日6：00~18：00之间每隔2h的室外温度（℃）如实验表1所视。

实验表1 室外温度观测结果（℃）

时间h 室内温度t1 室外温度t2 6 18.0 15.0 8 20.0 19.0 10 22.0 24.0 12 25.0 28.0 14 30.0 34.0 16 28.0 32.0 18 24.0 30.0 试用三次样条插值分别求该日室内外6：00～17：30之间每隔2h各点的近似温度（℃）

流程图:

开始 x=[6 8 10 12 14 16 18]; y1=[18.0 20.0 22.0 25.0 30.0 28.0 24.0]; y2=[15.0 19.0 24.0 28.0 34.0 32.0 30.0]; x0=[6.30:2.0:17.30] 使用函数spline(x,y,x0)) 结束

源程序代码：

%三次样条插值函数：spline(x,y,x0) %整体时间：x %所抽取的时间：x0 %温度：y

x=[6 8 10 12 14 16 18]; %整体时间 y1=[18.0 20.0 22.0 25.0 30.0 28.0 24.0]; %室内温度 y2=[15.0 19.0 24.0 28.0 34.0 32.0 30.0]; %室外温度

x0=[6.30:2.0:17.30]; %外6：00～17：30之间每隔2h的时间 y11=spline(x,y1,x0) %使用三次样条插值函数 y12=spline(x,y2,x0) %使用三次样条插值函数

结果

3．已知lgx在[1，101]区间10个整数采样点的函数值如实验表2所示。

试求lgx的5次拟合多项式p(x)，并绘制出lgx和p(x)在[1，101]区间的函数曲线。

流程图：

开始 x=[1 11 21 31 41 51 61 71 81 91 101]; lgx=[0 1.0414 1.322 1.4914 1.6128 1.7076 1.7853 1.8513 1.9085 1.9590 2.0043]; 使用函数p=polyfit(x,lgx,5) yi=polyval(p,[0:5:100]); plot(x,lgx,'r*'); plot(xi,yi); 结束

源程序代码：

x=[1 11 21 31 41 51 61 71 81 91 101];

lgx=[0 1.0414 1.322 1.4914 1.6128 1.7076 1.7853 1.8513 1.9085 1.9590 2.0043]; p=polyfit(x,lgx,5) xi=[0:5:100];

yi=polyval(p,[0:5:100]); plot(x,lgx,'r*'); hold on; plot(xi,yi);

legend('采样点','polyfit拟合函数曲线'); xlabel('x'); ylabel('lgx'); title('函数曲线');

结果：

因此lgx的5次拟合多项式p(x)= -0.1326x^5+0.1536x^4-0.0056x^3-0.0001x^2 函数曲线：