高振荡函数的数值积分Levin方法【MATLAB代码】

function Levin % Compute the integration of highly oscilatory function % \int_0^1 cos(log(x)/x)/x dx % using Levin Method warning off %#ok a = 1e-32; % less than eps(1e-16) b = 1; n = 21; % # of nodes f = @(x)exp(x); %1./x; %% Chebyshev Tn(x) w = @(x) 80*x;%log(x)./x; dw = @(x) 80*ones(size(x));%(1 - log(x))./x.^2; fprintf('================================================\n') fprintf('基函数选择：第一类切比雪夫函数T_k(x)\n') fprintf('------------------------------------------------\n') fprintf('配置点个数: n = %d\n\n',n); fprintf('配置点位置 积分结果\n') tic; for iter=1:4 if iter==1 x=linspace(a,b,n).'; elseif iter==2 x = cos((n-1:-1:0).'*pi/n); % extreme points of Chebyshev polynomial x = scale(x,[min(x),max(x)],[a,b]); elseif iter==3 x = cos((2*(n-1:-1:0).'+1)/n*pi/2); % zero points of Chebyshev polynomial x = scale(x,[min(x),max(x)],[a,b]); % linearly sacle x from [-1 1] to [0 1] elseif iter==4 % pnodes = [ ... % 0.2077849550078985; 0.4058451513773972; 0.5860872354676911; ... % 0.7415311855993944; 0.8648644233597691; 0.9491079123427585; ... % 0.9914553711208126]; % x = [-pnodes(end:-1:1); 0; pnodes]; % x=[0.0021714184870960 0.0130467357414141 0.0349212543221459 0.0674683166555077 0.1095911367067916 0.1602952158504878 0.2186214326656977 0.2833023029353764 0.3528035686492699 0.4255628305091844 0.5000000000000000 0.5744371694908156 0.6471964313507301 0.7166976970646236 0.7813785673343023 0.8397047841495122 0.8904088632932084 0.9325316833444923 0.9650787456778541 0.9869532642585859 0.9978285815129040].'; x = scale(x,[min(x),max(x)],[a,b]); end [k,x]=meshgrid(0:n-1,x); u = chebyshevT(k,x); du = [zeros(n,1) k(:,2:end).*chebyshevU(k(:,2:end)-1,x(:,2:end))]; A = du+1i*dw(x).*u; alpha = A\f(x(:,1)); % coefficients I = real(chebyshevT(k(1,:),b)*alpha*exp(1i*w(b)) - chebyshevT(k(1,:),a)*alpha*exp(1i*w(a))); elapse=toc; if iter==1 fprintf('[0,1]上均匀分布的点 %.15g\n',I) elseif iter==2 fprintf('[0,1]上切比雪夫极值点 %.15g\n',I) elseif iter==3 fprintf('[0,1]上切比雪夫零点 %.15g\n',I) elseif iter==4 fprintf('[0,1]上Gauss-Kronrod点 %.15g\n',I) end end fprintf('\n总共耗时：%g secs.\n',elapse) fprintf('------------------------------------------------\n\n\n') %% power function % original code: winner245 q = w; qpm = dw; fprintf('================================================\n') fprintf('基函数选择：幂函数 (x-(a+b)/2)^k\n') fprintf('------------------------------------------------\n') fprintf('配置点个数: n = %d\n\n',n); fprintf('配置点位置 积分结果\n') tic; d = (a+b)/2+eps((a+b)/2); for iter=1:4 if iter==1 x=linspace(a,b,n).'; elseif iter==2 x = cos((n-1:-1:0).'*pi/n); % extreme points of Chebyshev polynomial x = scale(x,[min(x),max(x)],[a,b]); elseif iter==3 x = cos((2*(n-1:-1:0).'+1)/n*pi/2); % zero points of Chebyshev polynomial x = scale(x,[min(x),max(x)],[a,b]); elseif iter==4 % pnodes = [ ... % 0.2077849550078985; 0.4058451513773972; 0.5860872354676911; ... % 0.7415311855993944; 0.8648644233597691; 0.9491079123427585; ... % 0.9914553711208126]; % x = [-pnodes(end:-1:1); 0; pnodes]; % x=[0.0021714184870960 0.0130467357414141 0.0349212543221459 0.0674683166555077 0.1095911367067916 0.1602952158504878 0.2186214326656977 0.2833023029353764 0.3528035686492699 0.4255628305091844 0.5000000000000000 0.5744371694908156 0.6471964313507301 0.7166976970646236 0.7813785673343023 0.8397047841495122 0.8904088632932084 0.9325316833444923 0.9650787456778541 0.9869532642585859 0.9978285815129040].'; x = scale(x,[min(x),max(x)],[a,b]); end u = bsxfun(@power, x-d, 0:n-1); uprime = bsxfun(@times, bsxfun(@power, x-d, -1:n-2), 0:n-1); qprime = repmat(qpm(x),1,n); A = uprime+1i*qprime.*u; alpha = A\f(x); I = real(u(end,:)*alpha*exp(1i*q(b)) - u(1,:)*alpha*exp(1i*q(a))); elapse=toc; if iter==1 fprintf('[0,1]上均匀分布的点 %.15g\n',I) elseif iter==2 fprintf('[0,1]上切比雪夫极值点 %.15g\n',I) elseif iter==3 fprintf('[0,1]上切比雪夫零点 %.15g\n',I) elseif iter==4 fprintf('[0,1]上Gauss-Kronrod点 %.15g\n',I) end end fprintf('\n总共耗时：%g secs.\n',elapse) fprintf('------------------------------------------------\n\n\n') end function y=chebyshevT(n,x) if isscalar(n) n=repmat(n,size(x)); elseif isscalar(x) x=repmat(x,size(n)); elseif ~isequal(size(n),size(x)) error('MATLAB:chebyshevT:WrongInputSize','Same size required for non-scalar inputs.') end sz=size(x); Nmax=max(n(:)); T=zeros(numel(x),Nmax+1); T(:,1)=1; T(:,2)=x(:); for k=2:Nmax+1 %0 to Nmax T(:,k+1)=2*x(:).*T(:,k)-T(:,k-1); end idx=sub2ind(size(T),1:numel(x),n(:)'+1); y=reshape(T(idx),sz); end function y=chebyshevU(n,x) if isscalar(n) n=repmat(n,size(x)); elseif isscalar(x) x=repmat(x,size(n)); elseif ~isequal(size(n),size(x)) error('MATLAB:chebyshevU:WrongInputSize','The same size required for non-scalar inputs.') end if any(n(:)<0)||any(floor(n(:))~=n(:)) error('MATLAB:chebyshevU:WrongValue','n must be non-negative integers.') end sz=size(x); Nmax=max(n(:)); U=zeros(numel(x),Nmax+1); U(:,1)=1; U(:,2)=2*x(:); for k=2:Nmax+1 %1 to Nmax U(:,k+1)=2*x(:).*U(:,k)-U(:,k-1); end idx=sub2ind(size(U),1:numel(x),n(:)'+1); y=reshape(U(idx),sz); end function y = scale(x,A,B) y=(x*diff(B)-A(1)*B(2)+A(2)*B(1))/diff(A); end