梯形法数值积分

本文格式为Word版，下载可任意编辑——梯形法数值积分

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

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

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

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

1.梯形法数值积分

A．算法说明：

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

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

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

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

在MATLAB中编程实现的复合梯形公式的函数为：CombineTraprl.功能：复合梯形公式求函数的数值积分。调用格式：[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)>epsI2否否n=n+1;h=(b-a)/n;I1=I2;I2=0;fori=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.0001endn=1;h=(b-a)/2;I1=0;

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

fori=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));endendI=I2;step=n;

D．应用举例:

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

开始?421dx2x?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.29392x?12．矩形法数值积分的源程序

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

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

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

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

fori=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);endendI=I2;step=n;

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

?21x2dx

a.流程图

开始f=x^2a=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=[];fori=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];endfori=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)];abcda2b2

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

结果：

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

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

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

（4）zcj=

400397396389389386384384382379378377377376375374374373371371368367366365365364364364363362362

362361361360358357357355353352352352352351351351350349349348348348348347345345344343342342341341340339338338336335335335333330330330

329328328327325323322322317315315314314312311310310308303299287285277275267

xsxh=

85936815972618538995391483956

13637724344430512762829850597086107922648788432581632922041589112314374673674257846

47725100806648547527295523174549761717838211152336942881906174993365601969409635

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

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

时间h室内温度t1室外温度t2618.015.0820.019.01022.024.01225.028.01430.034.01628.032.01824.030.0试用三次样条插值分别求该日室内外6：00～17：30之间每隔2h各点的近似温度（℃）

流程图:

开始x=[681012141618];y1=[18.020.022.025.030.028.024.0];y2=[15.019.024.028.034.032.030.0];x0=[6.30:2.0:17.30]使用函数spline(x,y,x0))终止

源程序代码：

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

x=[681012141618];%整体时间y1=[18.020.022.025.030.028.024.0];%室内温度y2=[15.019.024.028.034.032.030.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=[1112131415161718191101];lgx=[01.04141.3221.49141.61281.70761.78531.85131.90851.95902.0043];使用函数p=polyfit(x,lgx,5)yi=polyval(p,[0:5:100]);plot(x,lgx,'r*');plot(xi,yi);终止

源程序代码：

x=[1112131415161718191101];

lgx=[01.04141.3221.49141.61281.70761.78531.85131.90851.9590