科學(xué)計算方法18(數(shù)值積分2)

1、插值觀點下的數(shù)值積分,多項式插值,分段低次多項式插值(分段線性插值和分段三次插值),Newton-Cotes家族,分段低階數(shù)值積分(復(fù)合數(shù)值積分),簡單函數(shù)(多項式、三角函數(shù)等)近似復(fù)雜函數(shù),基本思想類似于分段插值,將積分區(qū)間適當(dāng)剖分,在每個子區(qū)間上用低階積分公式,這種方法稱為復(fù)合(composite)積分公式。,復(fù)合梯形積分公式,n 等分積分區(qū)間a,b 記h=(b-a)/n, xj=a+jh (j=0,n),clear all f=inline(x.(0.5);%高階版本建議匿名函數(shù) a=0.1;%積分左端點 b=2;%積分右端點 n=10000; x=linspace(a,b,n);%n+

2、1等分積分區(qū)間 y=f(x); %左矩陣公式 %int_left=(b-a)/(n+1)*sum(y(1:end-1) %右矩陣公式 %int_right=(b-a)/(n+1)*sum(y(1:end-1) %梯形公式 int_trap=0.5*(b-a)/(n+1)*sum(y(1:end-1)+y(2:end),前面的做法,例：用梯形法計算下面定積分 ( 取 n=100 ),解：,a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 ),x=0:1/100:1; y=1./(1+x.2); inum=trapz(x, y),trapz函數(shù),trapz 舉例,f

3、=inline(x.(0.5);%高階版本建議匿名函數(shù) a=0.1;b=2;n=10000; x=linspace(a,b,n); y=f(x); figure,plot(x,y) n=4; x1=linspace(a,b,n); y1=f(x1); for i=1:length(x1)-1 hold on, fill(x1(i),x1(i+1),x1(i+1),x1(i),x1(i),0,0,y1(i+1),y1(i),0,r) end,基本思想類似于分段插值,將積分區(qū)間適當(dāng)剖分,在每個子區(qū)間上用低階積分公式,這種方法稱為復(fù)合(composite)積分公式。,復(fù)合梯形積分公式,一般介值定理,

4、設(shè)函數(shù)y=f(x)在閉區(qū)間a,b上連續(xù),x1,xn是閉區(qū)間a,b中的點。對a1,an0, 在a和b之間存在滿足(a1+an) f()= a1 f(x1) +an f(xn)。,記 h=(b a)/(2n) , xj=a+jh ( j =0,1,2n),復(fù)合Simpson公式,Simpsion公式,quad(f,a,b,tol) f = f(x) 為被積函數(shù), a,b 為積分區(qū)間, tol 為計算精度,不用自己分割積分區(qū)間 可以指定計算精度, 若不指定, 缺省精度是 10-6 精度越高, 函數(shù)運行的時間越長 f 是函數(shù)句柄, 其中涉及的運算必須采用向量化,quad,自適應(yīng)Simpson方法,解：

5、,f=(x) 1./(1+x.2); inum=quad(f, 0, 1) % 采用缺省精度,inum=quad(x) 1./(1+x.2), 0, 1, 1e-10),例：用 quad 計算定積分：,quad 舉例,integral(f,a,b) integral(f,a,b,RelTol,tol),該函數(shù)比 quad 效率更高, 且可以處理一些非正常積分 可以指定計算精度, 若不指定, 缺省精度是 10-6 f 必須是函數(shù)句柄, 且涉及的運算必須采用向量化,integral,全局自適應(yīng)積分法（R2012a以后版本）,f=(x) 1./(1+x.2); inum=integral(f,0,1

6、) inum=integral(f,0,1,RelTol,1e-10),f=(x) exp(-x); inum=integral(f,0,inf),integral2(f,a,b,c,d,tol) integral2(f,a,b,c,d,RelTol,tol),可以指定計算精度, 若不指定, 缺省精度是 10-6 f 必須是函數(shù)句柄, 且涉及的運算必須采用向量化,integral2,計算二重積分的全局自適應(yīng)積分法,凡線面體皆設(shè)為由小漸大,一剎那中所增之積即微分也。其全積即積分也,蒙特卡洛(Monte Carlo method)方法,N=1000; x=rand(1,N);y=rand(1,N)

7、; plot(x,y,.),Ulam: 一個數(shù)學(xué)家的經(jīng)歷,蒙特卡洛(Monte Carlo method)方法,N=100000; data=rand(N,2); x=data(:,1)+1; y=data(:,2)*sin(1); II=find(y=sin(x)./x); M=length(II); S=M/N*sin(1),MonteCarlo積分的最大優(yōu)勢就在于高維積分，以及不規(guī)則區(qū)域。,蒙特卡洛(Monte Carlo method)方法,N=1000; x=rand(1,N);y=rand(1,N);z=rand(1,N); plot3(x,y,z,.),fun=(x,y)1+sqrt(1-x.2-y.2)-sqrt(x.2+y.2) ymax=(x)sqrt(1-x.2); ymin=(x)-sqrt(1-x.2); integral2(fun,-1,1,ymin,ymax),x,y=meshgrid(-1:0.001:1,-1:0.001:1); z=1+sqrt(1-x.2-y.2)-sqrt(x.2+y.2); mesh(x,y,real(z) ind=find(x.2+y.2=1); max(z(ind) min(z(ind),練習(xí)：,Reference: ,x=unifrnd(-1,1,100000,1); y=unifrnd(-1,1,1000