数值积分及matlab实现综述.doc

《数值积分及matlab实现综述.doc》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。

1、数值积分及MATLAB实现综述各种求积公式的MATLAB编程实现与应用MATLAB是由MathWorks公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过MATLAB软件编程实现，以下程序均用MATLAB7.0编写，运行坏境：1、硬件环境CPU(intelCorei3-2310M,2.1GHz),内存(2GB昱联)，2、软件环境w

2、indows7(32位)操作系统。以下总共编写了六个算法程序，部分代码参考文献[10-13],为了体现程序的正确性，以下程序都以/=「竺为例进行运算。原积分的精确值为Jox■/=f^4仕=0.946083070367183牛顿-科特斯求积公式的MATLAB实现先用M文件定义一个名为fl.m的函数：%i是要调用第几个被积函数g(i),x是自变量functionf=fl(i,x)g(l)=sqrt(x);ifx=()g(2)=l；elseg(2)=sin(x)/x;endg⑶=4/(l+xA2);f=g(i)；程序一:function[C,g]二NCot

3、es(a,b,n,m)%a,b分别为积分的上下限；%n是子区间的个数；%m是调用上面第儿个被积函数；%当n二1时计算梯形公式；当n二2时计算辛浦牛公式，以此类推;i=n;h=(b-a)/i;z=0;forj=0:ix(j+l)二a+j*h;ifj==0elsefork二1:js=s*k;endendr=l;ifi-j--0r=r;elsefork=l:(i-j)r=r*k;endendifmod((i-j),2)==lq二(i*s*r);elseq=i*s*r;fork=0:iifk〜二jy=y*(sym(,t,)-k);endendl=int(y,

4、0,i);C(j+l)=l/q；z=z+C(j+1)*f1(m,x(j+1));endg=(b-a)*z1)当输入a=O,b=l,刃=1,加=2时，即在MATLAB命令窗口输入»NCotes(0,l,l,2)即可得用梯形公式的积分值和相应科特斯系数如图3.11)当输入a=O,b=,巾=2,加=2时，即在MATLAB命令窗口输入»NCotes(0,l,2,2)即可得用辛浦生公式的积分值和相应科特斯系数如图3.22)当输入a=O,b=,斤=4,加=2时，即在MATLAB命令窗口输入»NCotes(0,l,4,2)即可得用科特斯公式的积分值和相应科特斯

5、系数如图3.3CommandWindow»NCotes(0,1,1,2)g=8293248040994423/9007199254740992arts=[1/2,1/2]»8293248040994423/9007199254740992arts=0.92073549240395»1卜—-图3.1图3.2图3.3复化求积公式的MATLAB实现一、复化梯形求积公式的MATLAB实现通过/（兀）的乳+1个等步长节点逼近积分嗨（加+他）+丐他）其中，xk=a^kh,x()=a,xn=ho程序二functions=trapr1(f,a,b,n)%f是被积函数

6、；%a,b分别为积分的上下限；%n是子区间的个数；%s是梯形总面积；h=(b-a)/n;s=0;fork=l:(n-l)x二a+h*k;s=s+feval('f,x);endformatlongs=h*(feval(,f,a)+feval('f,b))/2+h*s;先用M文件定义一个名为f.m的函数:functiony=f(x)ifx=0y二1;elsey=sin(x)/x;end在MATLAB命令窗口中输入»traprl(T,0,l,4)回车得到如图3.4图3.4若取子区间的个数/1=8在MATLAB命令窗口中输入»traprl(T,0,l,8)回

7、车得到如图3.5图3.5龙贝格求积公式的MATLAB实现构造T数表来逼近积分[fx)dx=R(J,J)其中。R(J,J)表示7数表的最后一行，最后一列的值。程序五：function[R,quad,err,h]=romber(f,a,b,n,delta)%f是被积函数%a,b分别是积分的上下限%n+1是T数表的列数%delta是允许误差%R是T数表%quad是所求积分值M=l;h=b-a;err=lJ二0;R=zeros(4,4)；R(1J)=h*(fevalCf,a)+feval(T,b))/2while((err>delta)&(J

8、(J

9、<4)J二J+l;h二h/2;s=();forp=l:Mx=a+h*(2*p-l);s=s+feval('f