梯形法数值积分

《MATLAB程序设计实践》课程考核一、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“矩形法、梯形法数值积分”1.梯形法数值积分A．算法说明：梯形法数值积分采用的梯形公式是最简单的数值积分公式，函数()fx在区间[a,b]上计算梯形法数值积分表达式为：()[()()]2babafxdxfafb由于用梯形公式来求积分十分粗糙，误差也比较大，后来改进后提出了复合梯形公式：bahn,其中，n为积分区间划分的个数；h为积分步长。在MATLAB中编程实现的复合梯形公式的函数为：CombineTraprl.功能：复合梯形公式求函数的数值积分。调用格式：[I,step]=CombineTraprl(f,a,b,eps).其中，f为函数名；a为积分下限；b为积分上限；eps为积分精度；I为积分值；Step为积分划分的区间个数B．流程图开始abs(I2–I1)epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1x=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;结束是否否n=1;h=(b-a)/2;I1=0;I2=(subs(sym(if),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;是C．复合梯形公式的MATLAB代码：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)/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．应用举例:复合梯形法求数值积分应用举例，利用复合梯形法计算定积分dxx42211a.流程图b.原程序代码：[q,s]=CombineTraprl('1/(x^2-1)',2,4)%精度为默认的10-4结果调用函数CombineTraprl(f,a,b,eps)f=1/(x^2-1)a=2;b=4结束开始[q,s]=CombineTraprl('1/(x^2-1)',2,4,1.0e-6)%精度为10-6结果：所以从复合梯形公式可以得出dxx422110.29392．矩形法数值积分的源程序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;应用举列：复合矩形法求数值积分应用举例，利用复合矩形法计算定积分dxx212a.流程图b.原程序代码：[q,s]=CombineTraprl2('x^2',1,2)%精度为默认的10-4结果：调用函数CombineTraprl2(f,a,b,eps)f=x^2a=1;b=2结束开始二、科学计算和工程实际问题和举例1.（题目）将100个学生5门功课的成绩存入矩阵P中，进行如下处理：（1）分别求每门课的最高分、最低分及相应学生序号。（2）分别求每门课的平均分和标准方差。（3）5门课总分的最高分、最低分及相应学生序号。（4）将5门课总分按从大到小顺序存入zcj中，相应学生序号存入xsxh。流程图：源程序代码：clearnum=[1:1:100]';p=[num,floor(50*rand(100,5)+45)];a=[];b=[];fori=1:5m=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:5m=min(p(:,i+1));pp=find(p(:,i+1)==min(p(:,i+1)));bb=ones(size(pp))*m;b=[b;0,0;pp,bb];使用函数max求出每门课最大值，并用find函数找到，然后输入到a；同理用min求出每门课最小值。使用mean函数求出每门课平均分c使用var函数求出每门课的标准方差使用sum函数求和，并找到每门课的最高分放入a2矩阵,每门课最低分放入b2。最后用sort函数和descend把总分排序放人zcj文件和xsxh文件构造一个矩阵做学生序号再构造一个随机矩阵做学生成绩让每门课最高分放入a矩阵，最低分放人b矩阵开始结束endc=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=400397396389389386384384382379378377377376375374374373371371368367366365365364364364363362362362361361360358357357355353352352352352351351351350349349348348348348347345345344343342342341341340339338338336335335335333330330330329328328327325323322322317315315314314312311310310308303299287285277275267xsxh=8593681597261853899539148395613637724344430512762829850597086107922648788432581632922041589112314374673674257846477251008066485475272955231745497617178382111523369428819061749933656019694096352．某气象观测站测得某日6：00~18：00之间每隔2h的室外温度（℃）如实验表1所视。实验表1室外温度观测结果（℃）时间h681012141618室内温度t118.020.022.025.030.028.024.0室外温度t215.019.024.028.034.032.030.0试用三次样条插值分别求该日室内外6：00～17：30之间每隔2h各点的近似温度（℃）流程图:源程序代码：%三次样条插值函数：spline(x,y,x0)%整体时间：x%所抽取的时间：x0%温度：yx=[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)%使用三次样条插值函数使用函数spline(x,y,x0))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]结束开始结果3．已知lgx在[1，101]区间10个整数采样点的函数值如实验表2所示。试求lgx的5次拟合多项式p(x)，并绘制出lgx和p(x)在[1，101]区间的函数曲线。流程图：使用函数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.95902.0043];开始结束源程序代码：x=[1112131415161718191101];lgx=[01.04141.3221.49141.61281.70761.78531.85131.90851.95902.0043];p=polyfit(x,lgx,5)xi=[0:5:100];yi=polyval(p,[0:5:100]);plot(x,lgx,'r*');holdon;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函数曲线：