数值计算&Matlab——高斯型积分

function varargout=gausscoef(a,b,n)% 高斯求积公式系数求解if nargin<3  n=2;    % 默认高斯点end x=sym('x%d',[1,n],'real');    % 节点A=sym('A%d',[1,n],'real');    % 系数eq=sym(zeros(2*n,1));         % 存储方程组for k=1:2*n  eq(k,1)=dot(A,x.^(k-1)-(b^k-a^k)/k);    % 构造方程endvars=symvar(eq);          % 提取eq的符号变量fun=matlabFunction(eq,'vars',{vars});   % 将eq转换成匿名函数的形式options=optimoptions(@fsolve,'FunctionTolerance',1e-12);x0=fsolve(fun,2*rand(1,2*n)-1,options);     % 求Gauss点和Gauss系数Ac=x0(1:n);     % Gauss系数x0(1:n)=[];     % Gauss点[varargout{1:2}]=deal(x0,...            % 第一个输出参数为Gauss点      Ac);                              % 第二个输出参数为Gauss求积系数

a = 2;b = 5;n=2;[x A]=gausscoef(a,b,n)

function p=opfun(polytype,n,xi) % opfun常用正交多项式% 存储各阶正交多项p=cell(n+1,1); % 递推多项式第二项的系数    P1=[1,0;1,0;2,0;-1,1;2,0]; % 判断所要求的多项式类型   ind=strcmpi(polytype,{'P','T','U','L','H'});% 递推多项式的首项p{1}=1;% 递推多项式的第二项p{2}=P1(ind,:);for k=2:n  switch polytype    case {'P','p'}   % 勒让德多项式      p{k+1}=(2*k-1)/k*[p{k},0]-(k-1)/k*[0,0,p{k-1}]; % 递推公式    case {'T','t','U','u'}     % 切比雪夫多项式      p{k+1}=2*[p{k},0]-[0,0,p{k-1}];    % 递推公式    case {'L','l'}             % 拉盖尔多项式      p{k+1}=conv([-1,2*k-1],p{k})-(k-1)^2*[0,0,p{k-1}];   % 递推公式    case {'H','h'}             % 埃尔米特多项式      p{k+1}=2*[p{k},0]-2*(k+1)*[0,0,p{k-1}];   % 递推公式      otherwise        % 如果多项式类型不是以上5种情形，则给出错误提示        error('无效的正交多项式类型.')      end     end if nargin==3               % 如果输入参数有3个  for k=1:n+1    p_val(k,:)=polyval(p{k},xi);   % 根据求得的正交多项式计算指定点处的值  end   p=p_val;                 % 将p_val赋给变量pend

function varargout = Gauss_1(n,type)% 常用高斯求积公式的高斯点及求积系数if nargin<2    type = 'P';    % 默认求Gauss-Legendre求积节点和求积系数endp=opfun(type(1),n);  % 由type指定正交多项式的系数p=p{end};            % 正交多项式P_n+1(x)x=roots(p);           % 求正交多项式的零点，即Gauss点switch type    case {'P',p}     % Gauss-Legendre求积        A=2./((1-x.^2).*polyval(polyder(p),x).^2);  % 求求积系数    case {'T','t'}   % Gauss-Chebyshev求积        A=pi/n*ones(size(x));                       % 求求积系数    case {'L','l'}   % Gauss-Laguerre求积        A=gamma(n+1)^2./(x.*polyval(polyder(p),x).^2);  % 求求积系数    case {'H','h'}    % Gauss-Hermite求积        A=2^(n+1)*gamma(n+1)*sqrt(pi)./polyval(polyder(p),x).^2;    case {'Pl','pl'}  % Gauss-Lobatto求积        x=roots(polyder(p)); % 求P'(x)的零点        A=[2/(n*(n+1));2./(n*(n+1)*polyval(p,x).^2);2./(n*(n+1))];        x=[-1;x;1];    % Gauss点    otherwise        error('本函数不能求其他类型的Gauss点及求积系数')end[varargout{1:2}]=deal(x,...                % 第一个输出参数为高斯点    A);                                    % 第二个输出参数为求积系数

f1 = @(x)x.^2./sqrt(1-x.^2);f2 = @(x)x.^2;[x01,A1]=Gauss_1(3,'P');[x02,A2]=Gauss_1(3,'T');I1=dot(A1,f1(x01))   % 高斯-勒让德求积I2=dot(A2,f2(x02))   % 高斯-切比雪夫求积

a=1;b=5;[x0,A]=Gauss_1(4,'Pl');f3=@(x) sin(x)./x;I3=(b-a)/2*dot(A,f3(x0*(b-a)/2+(a+b)/2))  % Gauss- Lobatto求积