《数值积分及MATLAB实现综述》由会员分享，可在线阅读，更多相关《数值积分及MATLAB实现综述（10页珍藏版）》请在金锄头文库上搜索。
1、 数值积分及MATLAB实现综述 各种求积公式的MATLAB编程实现与应用MATLAB是由MathWorks公式开发的一种主要用于数值计算及可视化图形处理的工程语言，是当今最优秀的科技应用软件之一。它将数值计算、矩阵运算、图形图像处理、信号处理和仿真等诸多强大的功能集成在较易使用的交互计算机环境中，为科学研究、工程应用提供了一种功能强、效率高的编程工具。下面我们将各种求积算法通过MATLAB软件编程实现，以下程序均用MATLAB7.0编写，运行坏境：1、硬件环境CPU（intel Core i3-2310M,2.1GHz），内存（2GB昱联）,2、软件环境windows7（32位）操作系统。以
2、下总共编写了六个算法程序，部分代码参考文献10-13,为了体现程序的正确性，以下程序都以为例进行运算。原积分的精确值为牛顿-科特斯求积公式的MATLAB实现先用M文件定义一个名为f1.m的函数：% i是要调用第几个被积函数g(i)，x是自变量function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x2);f=g(i);程序一：function C,g=NCotes(a,b,n,m)% a，b分别为积分的上下限；% n是子区间的个数；% m是调用上面第几个被积函数；% 当n=1时计算梯形公式；当n=
3、2时计算辛浦生公式，以此类推； i=n; h=(b-a)/i; z=0;for j=0:i x(j+1)=a+j*h; s=1; if j=0 s=s; elsefor k=1:j s=s*k;endendr=1;if i-j=0 r=r;elsefor k=1:(i-j) r=r*k;endendif mod(i-j),2)=1 q=-(i*s*r);else q=i*s*r;endy=1;for k=0:i if k=j y=y*(sym(t)-k); endendl=int(y,0,i);C(j+1)=l/q; z=z+C(j+1)*f1(m,x(j+1);endg=(b-a)*z1）当
4、输入，时，即在MATLAB命令窗口输入 NCotes(0,1,1,2)即可得用梯形公式的积分值和相应科特斯系数 如图3.12）当输入，时，即在MATLAB命令窗口输入 NCotes(0,1,2,2)即可得用辛浦生公式的积分值和相应科特斯系数如图3.23）当输入，时，即在MATLAB命令窗口输入 NCotes(0,1,4,2)即可得用科特斯公式的积分值和相应科特斯系数如图3.3图 3.1 图 3.2 图3.3 复化求积公式的MATLAB实现一、复化梯形求积公式的MATLAB实现通过的个等步长节点逼近积分其中，。程序二：function s=trapr1(f,a,b,n)% f是被积函数；% a,
5、b分别为积分的上下限；% n是子区间的个数；% s是梯形总面积；h=(b-a)/n;s=0;for k=1:(n-1) x=a+h*k; s=s+feval(f,x);endformat long s=h*(feval(f,a)+feval(f,b)/2+h*s;先用M文件定义一个名为f.m的函数：function y=f(x)if x=0 y=1;else y=sin(x)/x;end在MATLAB命令窗口中输入 trapr1(f,0,1,4) 回车得到 如图3.4 图 3.4若取子区间的个数在MATLAB命令窗口中输入 trapr1(f,0,1,8) 回车得到 如图3.5 图3.5 龙贝格
6、求积公式的MATLAB实现构造数表来逼近积分其中。表示数表的最后一行，最后一列的值。程序五：function R,quad,err,h=romber(f,a,b,n,delta)% f是被积函数% a,b分别是积分的上下限% n+1是T数表的列数% delta是允许误差% R是T数表% quad是所求积分值M=1;h=b-a;err=1J=0;R=zeros(4,4);R(1,1)=h*(feval(f,a)+feval(f,b)/2while (errdelta)&(Jn)|(J romber(f,0,1,5,0.5*(10(-8)回车得到 如图3.6 图3.6 高斯-勒让德求积公式的MAT
7、LAB实现程序六：function A,x=Guass1(N) i=N+1; f=(sym(t)2-1)i; f=diff(f,i); t=solve(f); for j=1:i for k=1:i X(j,k)=t(k)(j-1); end if mod(j,2)=0 B(j)=0; else B(j)=2/j; end end X=inv(X); for j=1:i A(j)=0; x(j)=0; for k=1:i A(j)=A(j)+X(j,k)*B(k); x(j)=x(j)+t(j); end x(j)=x(j)/k; endfunction g= GuassLegendre (a
8、,b,n,m)% a,b分别是积分的上下限；% n+1为节点个数；% m是调用f1.m中第几个被积函数； A,x=Guass1(n);g=0;for i=1:n+1 y(i)=(b-a)/2*x(i)+(a+b)/2; f(i)=f1(m,y(i); g=g+(b-a)/2*f(i)*A(i);end用M文件分别把上面两个自定义函数定义为名为Guass1.m函数和GuassLegendre.m函数用M文件定义一个名为f1.m的函数function f=f1(i,x) g(1)=sqrt(x);if x=0 g(2)=1;elseg(2)=sin(x)/x;endg(3)=4/(1+x2);f=g(i);在MATLAB命令窗口中输入 GuassLegendre (0,1,2,2) GuassLegendre (0,1,3,2)回车得到 如图3.7 图3.7- 9 -
