使用 Gauss-Legendre 求积进行双重积分的更快代码

计算科学 matlab 数值分析 正交
2021-12-07 04:02:45

我提供了以下代码来评估在 MatLab 中使用高斯勒让德求积的双积分

m=100;

%generate weights and abscissas
[wx,xx]=leg(-1,1,m);
[wtheta,xtheta]=leg(0,2.*pi,m);

%define function
psi=@(x,theta) hypergeom(-3./4,1./2,x.^2.*exp(1i.*theta));

%integrate with respect to x
intx=zeros(1,m);
for num=1:m
    intx(num)=sum(wx.*psi(xx,xtheta(num)));
end

sum(wtheta.*intx)

我在不同的脚本中定义了函数 leg(x1,x2,m) 来生成权重和横坐标，我只是在我的代码中调用它。与 Mathematica 的 NIntegrate 相比，我的 MatLab 代码运行缓慢。

我想让我的代码更快，因为我使用 MatLab 的想法是它比 Mathematica 更快。有什么办法可以让我的代码运行得更快？

附件是我运行代码时的个人资料摘要

资料概要

1个回答

原因可能是您的方法没有那么复杂。我敢打赌 NIntegrate 使用自适应方法并将积分分成几部分，然后根据误差估计进行细化。如果你不改进你的算法，你可能不会打败它。正如有人指出的那样，mupadex 是昂贵的部分。这来自函数评估的数量（特别是 hypergeom）。您可以通过一些分而治之的算法来减少它们。无论如何，使用如此多的高斯节点并不是一个好主意。更好的方法是使用低阶积分器并细化网格。如果您需要代码，请查看 matlab 函数integral2Calc.m 以了解实际算法（integral2t）或其包装器integral2.m

如果有人想运行上面的代码，这里有一个工作腿函数。

function [x, w] = leg(a,b,N)
% Generates the abscissa and weights for a Gauss-Legendre quadrature.
beta = (1:N-1)./sqrt(4*((1:N-1)).^2-1);
[w,x] = eig(diag(beta,-1)+diag(beta,1));
absc = diag(x);
wght = 2*w(1,:)'.^2;

% Linear map from[-1,1] to [a,b]
x=(a*(1-absc)+b*(1+absc))/2;      
w=wght;

x=x';
w=w';

end
其它你可能感兴趣的问题
上一篇FEM接触问题的新手帮助 下一篇基于测量序列生成伪随机序列