数值积分matlab程序.doc

第二章 数值积分2.1. 复化 Simpson 公式功能：利用复化 Simpson 公式计算被积函数 f(x)在给定区间上的积分值-function S=FSimpson(f,a,b,n)% f:被积函数句柄% a,b:积分区间的两个端点% n:子区间个数% S:用复化 Simpson 法求得的积分值h=(b-a)/n;fa=feval(f,a);fb=feval(f,b);S=fa+fb;x=a;for i=1:Nx=x+h/2;fx=feval(f,x);S=S+4*fx;x=x+h/2;fx=feval(f,x);S=S+2*fx;endS=h*S/6;-附：函数值为向量形式的 simpson 求积法function I=simpson_h(f,h)%调用格式 Isimpson(f,h)%f 为一向量，指定已知节点处的函数值%h 为步长n=length(f)-1;if n=1fprintf(Data has only one interval),return;end;if n=2I=(h/3)*(f(1)+4*f(2)+f(3);return;end;if n=3I=(3*h/8)*(f(1)+3*f(2)+3*f(3)+f(4);return;end;I=0;if 2*floor(n)=n % floor is a function round towards -infI=3*(h/8)*(f(n-2)+3*f(n-1)+3*f(n)+f(n+1);m=n-3;elsem=n;end;I=I+(h/3)*(f(1)+4*sum(f(2:2:m)+f(m+1);if m2I=I+(h/3)*2*sum(f(3:2:m);end;-附：函数值为向量形式的复合 simpson 求积法function I=simpson_n(fname,a,b,n)%调用格式： I=simpson_n(fname,a,b,n)%其中 a,b 为积分区间两个端点， n 为子区间数目h=(b-a)/n;x=a+(0:n)*h;f=feval(fname,x);I=simpson_h(f,h) % 调用上面编译好的 simpson_h 函数-2.2. 变步长梯形法功能：利用变步长梯形法计算函数 f(x)在给定区间的积分值-function T,n=bbct(f,a,b,eps)% f:被积函数句柄% a,b:积分区间的两个端点% eps:精度% n:二分区间的次数% T:用变步长梯形法求得的积分值h=b-a;fa=feval(f,a);fb=feval(f,b);T1=h*(fa+fb)/2;T2=T1/2+h*feval(f,a+h/2)/2;n=1;%按变步长梯形法求积分值while abs(T2-T1)=epsh=h/2;T1=T2;S=0;x=x+h/2;while xepsJ=J+1;h=h/2;S=0;for p=1:Mx=a+h*(2*p-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;M=2*M;for k=1:JR(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k)/(4k-1);enderr=abs(R(J+1,J)-R(J+1,J+1);endquad=R(J+1,J+1);2.4. 三点 Gauss 公式功能：利用三点 Gauss 公式计算被积函数 f(x)在给定区间的积分值-function G=TGauss(f,a,b)% f:被积函数句柄% a,b:积分区间的两个端点% G:用三点 Gauss 公式法求得的积分值x1=(a+b)/2-sqrt(3/5)*(b-a)/2;x2=(a+b)/2+sq