高震荡函数的积分

winner245

先说明一下，这个帖子的讨论，来源于MATLAB中文论坛：http://www.ilovematlab.cn/thread-264414-1-1.html

鉴于本论坛高手众多，只是过于冷清，所以，想把这个问题带到这里讨论一下。还请各位版主、高手不吝赐教。先说说问题描述。

如何求解下列数值积分？这个被积函数在0附近高度震荡，所以，数值求解也相当麻烦，直接用 integral 函数是求不出解的


不知道，这类积分该如何精确求解

winner245

我先抛砖引玉，给出我自己的解法：

1. format long
   
2. warning off
   
3. tol = 1e-6;  % 6 位小数精度
   
4. f = @(x) cos(log(x)./x)./x;
   
5. [y,e] = quadgk(f,1e-2,1,'RelTol',1e-8,'AbsTol',1e-12);
   
6. error = [];
   
7. ub = 1e-2; lb = ub/2;
   
8. [y0,e] = quadgk(f,lb,ub,'RelTol',1e-8,'AbsTol',1e-12);
   
9. while abs(y0)>tol
   
10.     if e>tol/100
    
11.         lb = (lb+ub)/2;
    
12.     else
    
13.         ub = lb;
    
14.         lb = ub/2;
    
15.         y = y+y0;    %　每次只叠加满足精度要求的积分
    
16.         error = [error e];  %记录一下每次有效积分的误差
    
17.     end  
    
18.     [y0,e] = quadgk(f,lb,ub,'RelTol',1e-8,'AbsTol',1e-12);
    
19. end
    
20. error
    
21. sum(error)
    
22. y

复制代码

得到：
y =
   0.323383328325528

liuyalong008

来一个下下策

1. 
   
2. numeric::quadrature(cos(log(x)/x)/x, x = 0..1, GaussLegendre = 10000)

复制代码

时间要很久

1. 0.3233675887181222334648423061006770509986438341236031749966615501914392850834872089215342986486185492

复制代码

winner245

liuyalong008 发表于 2013-10-15 15:14
来一个下下策时间要很久

呵呵，这个不错，关键是精度很高。这个是mupad吧？

bainhome

样条针对高震荡积分应该是个不错的策略，我电脑上没装MATLAB，谁有时间可以试试B样条，以前试过cos(17x)的导数，还是很精确的。

winner245

bainhome 发表于 2013-10-15 23:52
样条针对高震荡积分应该是个不错的策略，我电脑上没装MATLAB，谁有时间可以试试B样条，以前试过cos(17x)的 ...

您好，感谢您的建议。之前没做过样条积分，刚才看了 help 文档，按照你的建议写了一个代码，不知道是不是你说的：

1. format long
   
2. f = @(x) cos(log(x)./x)./x;
   
3. x = linspace(1e-10,1,1e5);
   
4. p = spline(x,f(x));
   
5. p = fn2fm(p,'B-');
   
6. diff(fnval(fnint(p),[0 1]))

复制代码

运行后结果为：
ans =
     3.343232137526005e+04

这个结果误差很大

bainhome

应该是spapi或csapi命令，总之一个是三次样条一个是B样条，这个貌似和spline算法不大一样。电脑不给力，不敢装MATLAB，无法测试想法，抱歉。

winner245

bainhome 发表于 2013-10-16 14:05
应该是spapi或csapi命令，总之一个是三次样条一个是B样条，这个貌似和spline算法不大一样。电脑不给力，不 ...

非常感谢！csapi 和 spline 都是3次样条，二者本质上一样。所以，我估计你说的是 spapi，这个是任意阶样条。我试了spapi 的20阶样条：

1. format long
   
2. f = @(x) cos(log(x)./x)./x;
   
3. x = linspace(1e-10,1,1e5);
   
4. k=20;
   
5. p = spapi(k,x,f(x));
   
6. diff(fnval(fnint(p),[0 1]))

复制代码

结果是：
ans =
     2.320685419330430e+04

也许我的使用方法不对...

winner245

liuyalong008 发表于 2013-10-15 15:14
来一个下下策时间要很久

请问你就是用上面原封不动的mupad代码算出了这么长的结果吗？我直接用你的代码计算不出来。另外，我试了mathematica：

1. NumberForm[NIntegrate[Cos[Log[x]/x]/x, {x, 0, 1}], 20]

复制代码

得到结果是：

0.3233674295759144

对比你的mupad结果

0.3233675887181222334648423061006770509986438341236031749966615501914392850834872089215342986486185492

差距很大，其中红色部分为不同的数字

liuyalong008

winner245 发表于 2013-11-2 03:30
请问你就是用上面原封不动的mupad代码算出了这么长的结果吗？我直接用你的代码计算不出来。另外，我试了ma ...

是用上述代码原封不动计算的，我的笔记本是没法计算的，用台式算的。至于差距较大，我个人还是觉得mathematica的结果应该更可靠，mupad如果要更高精度的话   时间要更久

winner245

liuyalong008 发表于 2013-11-4 12:12
是用上述代码原封不动计算的，我的笔记本是没法计算的，用台式算的。至于差距较大，我个人还是觉得mathem ...

嗯，最近发现mathematica在计算数值积分这一块很强大。对于震荡函数积分，mathematica里有全面的考虑，而且mathematica在考虑震荡积分时，就是拿我上面这个积分做示例去计算的。而maple里却找不到半点关于高震荡积分的资料。我最近一直在研究如何手动实现高震荡积分，等有结果了，我把matlab代码发上来。

winner245

这个问题终于知道怎么用matlab解了。

最近得空读了一些高震荡函数积分的算法，使用 Levin 算法可以解决问题。先附上论文：

代码如下

1. format long
   
2. warning off
   
3. a = 1e-20;
   
4. b = 1;
   
5. n = 20;
   
6. f = @(x) 1./x;
   
7. q = @(x) log(x)./x;
   
8. qpm = @(x) (1 - log(x))./x.^2;
   
9. d = (a+b)/2+eps((a+b)/2);
   
10. x = linspace(a,b,n).';
    
11. u = bsxfun(@power, x-d, 0:n-1);
    
12. uprime = bsxfun(@times, bsxfun(@power, x-d, -1:n-2), 0:n-1);
    
13. qprime = repmat(qpm(x),1,n);
    
14. A = uprime+1i*qprime.*u;
    
15. alpha = A\f(x);
    
16. I = real(u(end,:)*alpha*exp(1i*q(b)) - u(1,:)*alpha*exp(1i*q(a)))

复制代码

I =
   0.323367760463858

winner245

bainhome 发表于 2013-10-15 23:52
样条针对高震荡积分应该是个不错的策略，我电脑上没装MATLAB，谁有时间可以试试B样条，以前试过cos(17x)的 ...

最近读了些高震荡函数积分的论文，才发现样条内插是一种有效的求解高震荡积分的算法，最早的针对样条的算法是 Filon 算法，采用的是2阶样条，后来该算法被推广到了B样条

我在12L给出的解法是 Levin 算法，这种算法不是样条法，它本质是一种匹配算法。下面对这个算法做一个简要的定性介绍。

Levin 算法是先将高震荡积分的被积函数表示成另一个函数的微分，这样，高震荡积分可以马上通过莱布尼兹公式写出结果。所以，高震荡积分就转变为了解微分方程。但我们并不需要求出微分方程的所有解，而只需要一组特解即可。但特解也不好直接求。于是，需要用到近似。这里的近似是将特解表示成为 n 个线性无关的基函数的线性组合，n 个组合系数待定。将这组特解代入微分方程，并选取积分区间上n个特定的点，从而得到关于n个组合系数的n元线性方程组。线性方程组很容易求得。这个求解过程是让n个组合系数和取的n个点去匹配微分方程，所以，匹配法因此而得名。很显然，基函数和n个点的选取至关重要。比如，基函数必须满足线性无关性。在Levin 82年和96年两篇论文里，它用的都是不同次数的多项式。n个匹配点是等间隔选取的。在Levin后，出现了一些改进算法，比如，基函数选择为 Chebychev多项式

高震荡积分的水很深，Levin 的两篇paper发表于82和96年，而直到今天，mathematica最新版的高震荡积分算法依然是这两篇论文的算法。虽然，最近10年出现了一些新论文，但基本上都是小的改进，没有大的飞跃。

bainhome

winner245 发表于 2014-1-6 08:47
最近读了些高震荡函数积分的论文，才发现样条内插是一种有效的求解高震荡积分的算法，最早的针对样条的算 ...

一直没用MATLAB算，因为没想到太好的办法改进，直觉样条的逼近算法是在相对“小”的区域内局部构造合理逼近的内插点，积分计算实际仍然是面积累加，这道题目中的振荡区域实在是太小，小到利用内插样条点构造累加区域时，计算累积误差相较之下，“大”到难以忍受。
看了winner245先生对Livin算法的介绍，自被积函数入手转为求微分方程，因为一般情况下积分是弱形式，微分是强形式，很难在解积分问题往微分方程上去想，高明！
winner245解释通透明白，又佐以代码以例释疑，同样很赞！

szldh2005

非常感谢winner245介绍高震荡函数的积分方法，获益匪浅。

winner245

bainhome 发表于 2013-10-15 23:52
样条针对高震荡积分应该是个不错的策略，我电脑上没装MATLAB，谁有时间可以试试B样条，以前试过cos(17x)的 ...

最近又试了一下你先前说的样条积分，对于震荡核为 cos(w*x) 或更一般的 exp(j*w*x) ，在 w 为固定的很大的数时（比如 w = 1000），样条法确实可以精确计算。

1. format long
   
2. f = @(x) cos(1000*x);
   
3. x = linspace(0,1,1e5);
   
4. k=20;
   
5. p = spapi(k,x,f(x));
   
6. diff(fnval(fnint(p),[0 1]))  % 样条法结果
   
7. integral(f,0,1)  % integral函数计算结果

复制代码

ans =
     8.268795405347876e-04

ans =
     8.268795405311970e-04

这里样条法之所以能奏效，是因为震荡速度虽然很快，但再快也只是一个固定震荡速度 （w=1000），所以，只要将被积函数划分出足够多的样条区间，再用合适阶数的样条内插就能得到很高的近似精度。而在1L的积分里，震荡速度在0点为无穷，此时，无论划分多少个样条区间，都很难达到满意的精度。这也是为什么我在8L用了完全一样的样条积分算法却始终得不到满意的结果。

szldh2005

winner245:
您好！能否贴一下Levin96的论文？

winner245

szldh2005 发表于 2014-1-13 08:38
winner245:
您好！能否贴一下Levin96的论文？

szldh2005

非常感谢winner245!