數(shù)值積分matlab程序.doc

第二章 數(shù)值積分2.1. 復(fù)化 Simpson 公式功能：利用復(fù)化 Simpson 公式計(jì)算被積函數(shù) f(x)在給定區(qū)間上的積分值-function S=FSimpson(f,a,b,n)% f:被積函數(shù)句柄% a,b:積分區(qū)間的兩個(gè)端點(diǎn)% n:子區(qū)間個(gè)數(shù)% S:用復(fù)化 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;-附：函數(shù)值為向量形式的 simpson 求積法function I=simpson_h(f,h)%調(diào)用格式 Isimpson(f,h)%f 為一向量，指定已知節(jié)點(diǎn)處的函數(shù)值%h 為步長(zhǎng)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;-附：函數(shù)值為向量形式的復(fù)合 simpson 求積法function I=simpson_n(fname,a,b,n)%調(diào)用格式： I=simpson_n(fname,a,b,n)%其中 a,b 為積分區(qū)間兩個(gè)端點(diǎn)， n 為子區(qū)間數(shù)目h=(b-a)/n;x=a+(0:n)*h;f=feval(fname,x);I=simpson_h(f,h) % 調(diào)用上面編譯好的 simpson_h 函數(shù)-2.2. 變步長(zhǎng)梯形法功能：利用變步長(zhǎng)梯形法計(jì)算函數(shù) f(x)在給定區(qū)間的積分值-function T,n=bbct(f,a,b,eps)% f:被積函數(shù)句柄% a,b:積分區(qū)間的兩個(gè)端點(diǎn)% eps:精度% n:二分區(qū)間的次數(shù)% T:用變步長(zhǎng)梯形法求得的積分值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;%按變步長(zhǎng)梯形法求積分值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. 三點(diǎn) Gauss 公式功能：利用三點(diǎn) Gauss 公式計(jì)算被積函數(shù) f(x)在給定區(qū)間的積分值-function G=TGauss(f,a,b)% f:被積函數(shù)句柄% a,b:積分區(qū)間的兩個(gè)端點(diǎn)% G:用三點(diǎn) Gauss 公式法求得的積分值x1=(a+b)/2-sqrt(3/5)*(b-a)/2;x2=(a+b)/2+sq