困难振荡积分的数值积分方法

计算科学 matlab 正交 特殊功能
2021-12-09 20:49:30

我需要对下面的积分进行数值评估:

0sinc(xr)rE(r)dr

在哪里E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2),xR+λ,κ,ν>0. 这里K是第二类修正贝塞尔函数。在我的特殊情况下,我有λ=0.00313,κ=0.00825ν=0.33.

我正在使用 MATLAB,并且尝试了内置函数integralquadgk,这给了我很多错误(见下文)。我自然也尝试了许多其他的事情,例如按部分积分,以及对积分求和kxπ(k+1)xπ.

那么,您对接下来我应该尝试哪种方法有什么建议吗?

更新(添加的问题)
我阅读了@Pedro 链接到的论文,我认为这并不难理解。但是,我有几个问题:

  • 可以用吗xk作为基础元素ψk,在单变量莱文方法中描述?
  • 我可以改用 Filon 方法吗,因为振荡的频率是固定的?

示例代码
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89

ans =

3.3197e+06

4个回答

我编写了自己的积分器,quadcc它比具有奇异点的 Matlab 积分器处理得更好,并提供了更可靠的误差估计。

要将它用于您的问题,我执行了以下操作:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

该函数f现在是您的被积函数。请注意,我刚刚将任何旧值分配给x.

为了在无限域上积分,我应用了变量替换:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

g从0到1积分应该和f从0到积分一样. 不同的变换可能会产生不同的质量结果:数学上所有变换都应该给出相同的结果,但不同的变换可能会产生更平滑或更容易积分g的 s。

然后我调用我自己的积分器quadcc,它可以处理NaN两端的 s:

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

请注意,误差估计很大,即quadcc对结果没有太大信心。不过,看看这个函数,这并不奇怪,因为它的振荡值比实际积分高三个数量级。同样,使用不同的间隔变换可能会产生更好的结果。

您可能还想查看更具体的方法,例如this它涉及更多,但绝对是解决此类问题的正确方法。

正如 Pedro 所指出的,Levin 类型的方法是解决这类问题的最佳方法。

您可以访问 Mathematica 吗?对于这个问题,Mathematica 默认会检测并使用它们:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

下面是 x 值范围内的图:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

从 x = 0.5 到 x=10 绘图

您还可以手动指定要应用的特定 Levin 类型的方法,在这种情况下可以产生轻微的性能改进:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

有关Mathematica 中 Levin 类型方法的详细信息,请参阅文档

如果您无法访问 Mathematica,您可以按照 Pedro 的建议在 Matlab 中编写 Levin 类型(或其他专门的振荡)方法。

你使用 Matlab 的chebfun库吗?我刚刚了解到它在这里包含一个基本的 Levin 类型方法的实现该实现由 Olver(振荡正交领域的专家之一)编写。它不处理奇点、自适应细分等,但它可能正是您开始所需要的。

Pedro 推荐的转换是一个好主意。您是否尝试过使用 Matlab 的“quadgk”函数中的参数?例如,使用 Pedro 的变换,您可以执行以下操作:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
使用这个给我一个解:
-2184689.50220729
并且只需要 0.8 秒(使用上面提到的值:x=10)
Walter Gander 和 Walter Gautschi 有一篇关于使用 Matlab 进行自适应求积的论文您也可以使用代码(链接在这里