查看: 30506|回复: 56|关注: 0
打印 上一主题 下一主题

快乐扑克3淘宝: MATLAB高效编程之向量化积分

  [复制链接]

论坛优秀回答者

快乐扑克3如何跟对子 www.vqxik.tw 24

主题

1万

帖子

1638

最佳答案
  • 关注者: 943
跳转到指定楼层
1#
发表于 2013-10-22 03:30:23 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
本帖最后由 winner245 于 2013-10-24 01:54 编辑

前些天和吴兄(编著《MATLAB高效编程技巧与应用:25个案例分析》)闲谈中,谈到了如何提高数值积分的效率,达成的共识是利用向量化积分来实现高效率的数值积分。

众所周知,MATLAB为我们提供了两个非常易用的向量化积分函数:quadv 和 integral 函数。本帖的目的是探讨一下,如何利用这两个积分函数,解决三类不同的向量化积分问题。

我先抛砖引玉,在跟帖里给出我个人这些天的体会?;队蠹也斡胩致?,批评指正,希望能把这个话题继续下去?;队魏问祷值幕疤?/font>

论坛优秀回答者

24

主题

1万

帖子

1638

最佳答案
  • 关注者: 943
来自 2#
 楼主| 发表于 2013-10-22 03:38:07 | 只看该作者
本帖最后由 winner245 于 2015-7-17 12:26 编辑

关于单个函数向量化积分的探讨
winner245
向量化积分简介

给定一个函数f,我们可能需要计算一系列定积分。比如,i) 求解一个带参数函数f(x,t),x为自变量,t为参数,在M个不同参数取值t= t1, …, tM时在某固定区间[a,b]上的定积分;ii)对于同一个函数f(x),求解在N个不同区间上[a1,b1]、[a2,b2]、、[aN,bN]的定积分;iii)求解带参数函数f(x,t)M个不同参数取值和N个不同区间上的定积分。显然,无论是上述哪一种情形,我们都将面临求解一系列定积分。本文重点关心的是如何通过向量化运算,求解出该函数所有的定积分,且计算过程不涉及任何循环(以及循环的简写代码:arrayfun)。本文将能够实现这类积分功能的向量化运算称之为“向量化积分”。

本文假定所求积分不存在解析解,即符号积分失效。故将重点讨论如何通过数值积分(借助quad、quadl、quadgk、quadv、integral等函数)来实现向量化积分。

问题描述

在数值计算中,我们常?;崤龅较铝腥嗷治侍?。

问题1(带参数函数在固定区间上的定积分):求带参数函数f(x,t) 当参数tM个不同数值时在某固定区间[a,b]上的定积分:


                              
问题2(不带参数函数在多个区间上的定积分):求函数f(x)N个不同区间[an,bn], ,上的定积分:


问题3(带参数函数在多个区间上的定积分):求带参数函数f(x,t) 当参数tM个不同数值,在N个不同区间[an,bn], n = 1,2, …, N,上的定积分:


为了方便后面写代码,本文将考虑下面的被积函数:







其中,x为自变量,t为参数。

问题分析
要求解一系列积分,最容易想到的办法就是循环求解。比如,通过for循环实现每次循环求解一个定积分。另外,for循环代码还可以通过arrayfun函数来简写,这种“简写”本质上还是循环,只不过代码看起来更简洁,代码效率与for循环相当(依然很低)。但下文仍会给出基于arrayfun的代码,作为参照系来衡量向量化积分的优越性。

对于问题1(固定区间上的带参数函数的定积分),可以直接调用MATLAB向量化积分函数quadvintegral,这两个函数是MATLAB提供的仅有的两个支持向量化积分的函数。具体做法是:将被积函数f(x,t) 定义为可接受向量化参数t的关于x的函数句柄,然后直接调用quadvintegral,其中,integral函数要求将'ArrayValued' 参数设置为 true。

问题2、3所示的积分都是在变化的多个区间上求解。然而,无论是quadv还是integral函数,在求解向量化积分时,均要求积分区间是一个固定区间(即积分上、下限里不存在任何变量)。如果不对积分做必要的处理,是无法实现一次quadv或integral调用同时求解多个不同区间上的定积分的。本文的做法是:使用换元法,将多个不同区间上的积分转换到一个相同区间(如[0,1])上,即实现将积分上、下限里的变量转移到被积函数里,从而得到固定区间上带参数函数的定积分,这也就是问题1描述的一类积分。最后,直接调用quadv或integral函数即可。

因为quadvquad的向量化版本,我将会比较(向量化的)quadv和(非向量化的)quad的效率。同样的,integral'ArrayValued' 参数为 true时,是integral在不带'ArrayValued' 参数时的向量化版本,故我们会比较integral函数在'ArrayValued'参数为truefalse时的效率。本文不会去比较向量化的quadv和向量化的integral之间的效率,也不会比较非向量化的quadintegral的效率,因为此类比较的结果严重依赖于被积函数类型,比如,与被积函数是否震荡、是否有奇异点有关,即不同的函数可能得到不同的结果。感兴趣的朋友还可以把quad换成quadl、quadgk来参与相关比较。

本文的所有测试结果均基于“i7-3630QM + 8GB RAM”硬件平台和“Win8 64 bit + MATLAB 2012b 64 bit”软件平台。

问题1解决方法

问题1代码:求参数t0 ~ 10之间的随机数构成的2*6矩阵时(即M = 12),f(x, t)在区间[0, 1] 上对x的所有积分值
  1. t = 10*rand(2,6);
  2. f = @(x,t) (x.^2-1)./(x.^2+1).*cos(t.*x.*exp(-x));
  3. I_quad = arrayfun(@(t) quad(@(x) f(x,t), 0, 1), t)
  4. I_integral = arrayfun(@(t) integral (@(x) f(x,t), 0, 1), t)
  5. I_quadv = quadv(@(x) f(x,t), 0, 1)
  6. I_integralv = integral(@(x) f(x,t), 0, 1, 'ArrayValued', true)
复制代码

运行后得到:

I_quad =
  -0.2576    0.0450   -0.4421  -0.2076   -0.3398   -0.5510
  -0.5348   -0.5549   -0.4850  -0.5576   -0.5535    0.0471
I_integral =
  -0.2576    0.0450   -0.4421  -0.2076   -0.3398   -0.5510
  -0.5348   -0.5549   -0.4850  -0.5576   -0.5535    0.0471
I_quadv =
  -0.2576    0.0450   -0.4421  -0.2076   -0.3398   -0.5510
  -0.5348   -0.5549   -0.4850  -0.5576   -0.5535    0.0471
I_integralv =
  -0.2576    0.0450   -0.4421  -0.2076   -0.3398   -0.5510
  -0.5348   -0.5549   -0.4850  -0.5576   -0.5535    0.0471

可以看出,向量化和非向量化积分得到了相同的结果。下面对比一下代码效率。

问题1代码效率:M = 1e4(即M = 1e4个不同参数取值)
  1. t = 10*rand(1,1e4);
  2. f = @(x,t) (x.^2-1)./(x.^2+1).*cos(t.*x.*exp(-x));
  3. tic; I_quad = arrayfun(@(t) quad(@(x) f(x,t), 0, 1), t); t_quad = toc;
  4. tic; I_integral = arrayfun(@(t) integral(@(x) f(x,t), 0, 1), t); t_integral = toc;
  5. tic; I_quadv = quadv(@(x) f(x,t), 0, 1); t_quadv = toc;
  6. tic; I_integralv = integral(@(x) f(x,t), 0, 1, 'ArrayValued', true); t_integralv = toc;
  7. disp([' t_quad = ', num2str(t_quad),', t_quadv = ',num2str(t_quadv),', t_integral = ',num2str(t_integral) ,', t_integralv = ',num2str(t_integralv)])
  8. disp([' t_quad/t_quadv = ', num2str(t_quad/t_quadv)])
  9. disp([' t_integral/t_integralv = ', num2str(t_integral/t_integralv)])
复制代码
运行结果如下:

  t_quad = 4.5577, t_quadv = 0.018436,t_integral = 5.3036, t_integralv = 0.053569
  t_quad/t_quadv = 247.2156
  t_integral/t_integralv = 99.0046

可见,向量化的quadv比非向量化的quad效率提高了247倍,向量化的integral比非向量化的integral效率提高了99倍。

问题2,3解决方法

问题2,3均可以通过如下式子换元:


则积分(2)、(3)分别转化为:




            
以上两个积分都是在固定区间上的带参数积分。

问题2代码:假设N= 10, 即针对10个区间[0,1], [1, 2],[2, 3], …, [9, 10]分别求解定积分
  1. a = [0:9];  % 积分下限
  2. b = a + 1;  % 积分上限
  3. f = @(x) (x.^2-1)./(x.^2+1).*cos(x.*exp(-x));
  4. g = @(x) (b-a).*f(a+(b-a)*x);
  5. I_quad = arrayfun(@(a, b) quad(f, a, b), a, b)
  6. I_integral = arrayfun(@(a, b) integral (f, a, b), a, b)
  7. I_quadv = quadv(g, 0, 1)
  8. I_integralv = integral(g, 0, 1, 'ArrayValued', true)
复制代码

运行结果为:

I_quad=
   -0.5551 0.3387  0.7009  0.8414 0.9036  0.9352  0.9534 0.9649  0.9726  0.9780
I_integral=
   -0.5551 0.3387  0.7009  0.8414 0.9036  0.9352  0.9534 0.9649  0.9726  0.9780
I_quadv=
   -0.5551 0.3387  0.7009  0.8414 0.9036  0.9352  0.9534 0.9649  0.9726  0.9780
I_integralv=
   -0.5551 0.3387  0.7009  0.8414 0.9036  0.9352  0.9534 0.9649  0.9726  0.9780

可见,向量化代码得到了和循环代码相同的结果。表面上看起来,循环代码似乎更简洁(主要归因于arrayfun的特性),但向量化代码在效率上有显著优势。下面给出a=[0:1e4]的代码效率对比:

问题2代码效率:N= 1e4+1(即N= 1e4 + 1个定积分区间)
  1. a = [0:1e4];  % 积分下限
  2. b = a + 1;    % 积分上限
  3. f = @(x) (x.^2-1)./(x.^2+1).*cos(x.*exp(-x));
  4. g = @(x) (b-a).*f(a+(b-a)*x);
  5. tic; I_quad = arrayfun(@(a, b) quad(f, a, b), a, b); t_quad = toc;
  6. tic; I_integral = arrayfun(@(a, b) integral (f, a, b), a, b); t_integral = toc;
  7. tic; I_quadv = quadv(g, 0, 1); t_quadv = toc;
  8. tic; I_integralv = integral(g, 0, 1, 'ArrayValued', true); t_integralv = toc;
  9. disp([' t_quad = ', num2str(t_quad),', t_quadv = ',num2str(t_quadv),', t_integral = ',num2str(t_integral) ,', t_integralv = ',num2str(t_integralv)])
  10. disp([' t_quad/t_quadv = ', num2str(t_quad/t_quadv)])
  11. disp([' t_integral/t_integralv = ', num2str(t_integral/t_integralv)])
复制代码

运行结果为:
  t_quad = 1.4122, t_quadv =0.014844, t_integral = 4.7597, t_integralv = 0.060604
  t_quad/t_quadv = 95.1414
  t_integral/t_integralv = 78.5379

可见,向量化的quadv比非向量化的quad效率提高了95倍,向量化的integral比非向量化的integral效率提高了78倍。

问题3代码:假设M= 2,N= 6, 即针对2个不同的参数值tt=1, 5)和6个不同区间[0,1], [1, 2],[2, 3], …, [5, 6]求解定积分。
  1. a = [0:5];   % 积分下限
  2. b = a + 1;   % 积分上限
  3. t = [1 5].';
  4. f = @(x,t) (x.^2-1)./(x.^2+1).*cos(t*x.*exp(-x));
  5. g = @(x,t) bsxfun(@times, (x.^2-1)./(x.^2+1), cos(t*(x.*exp(-x))));
  6. g = @(x,t) bsxfun(@times, b-a, g(a+(b-a)*x, t));
  7. I_quad = cell2mat(arrayfun(@(t) arrayfun(@(a, b) quad(@(x) f(x,t), a, b), a, b), t, 'UniformOutput', false))
  8. I_integral = cell2mat(arrayfun(@(t) arrayfun(@(a, b) integral(@(x) f(x,t), a, b), a, b), t, 'UniformOutput', false))
  9. I_quadv = quadv(@(x) g(x,t), 0, 1)
  10. I_integralv = integral(@(x) g(x,t), 0, 1, 'ArrayValued', true)
复制代码

问题3代码效率:M = 100, N = 101(即M= 100个不同参数取值,N= 101个定积分区间)
  1. a = [0:100];   
  2. b = a + 1;   
  3. t = 10*rand(100,1);
  4. f = @(x,t) (x.^2-1)./(x.^2+1).*cos(t*x.*exp(-x));
  5. g = @(x,t) bsxfun(@times, (x.^2-1)./(x.^2+1), cos(t*(x.*exp(-x))));
  6. g = @(x,t) bsxfun(@times, b-a, g(a+(b-a)*x, t));
  7. tic; I_quad = cell2mat(arrayfun(@(t) arrayfun(@(a, b) quad(@(x) f(x,t), a, b), a, b), t, 'UniformOutput', false)); t_quad = toc
  8. tic; I_integral = cell2mat(arrayfun(@(t) arrayfun(@(a, b) integral(@(x) f(x,t), a, b), a, b), t, 'UniformOutput', false)); t_integral = toc
  9. tic; I_quadv = quadv(@(x) g(x,t), 0, 1); t_quadv = toc
  10. tic; I_integralv = integral(@(x) g(x,t), 0, 1, 'ArrayValued', true); t_integralv = toc
  11. disp([' t_quad = ', num2str(t_quad),', t_quadv = ',num2str(t_quadv),', t_integral = ',num2str(t_integral) ,', t_integralv = ',num2str(t_integralv)])
  12. disp([' t_quad/t_quadv = ', num2str(t_quad/t_quadv)])
  13. disp([' t_integral/t_integralv = ', num2str(t_integral/t_integralv)])
复制代码

运行结果为:

  t_quad = 2.2115, t_quadv =0.022581, t_integral = 5.4, t_integralv = 0.06808
  t_quad/t_quadv = 97.9367
  t_integral/t_integralv = 79.3189

可见,向量化的quadv比非向量化的quad效率提高了97倍,向量化的integral比非向量化的integral效率提高了79倍。

结论

本文重点介绍了向量化积分利器:quadv integral函数,以及如何利用这两个函数实现三类不同的向量化积分?;镜乃悸肥墙蠡肿?/font>quadvintegral函数可以接受的有效形式,然后直接利用这两个函数的求出一系列积分。通过比较,证实了在向量化尺寸足够大时(如1e4级别),向量化积分的效率远高于循环实现。

Future Topic:关于多个函数向量化积分的探讨

本文仅探讨了如何对单个函数(带参数f(x,t) 或者不带参数f(x))的向量化积分方法。本人计划在下一期分享帖里,继续探讨同时对多个不同函数f1,f2, f3, fn(带参数或不带参数)在任意个区间上的向量化积分方法,该方法将实现通过一次integralquadv调用计算出所有(带参数或不带参数)函数在任意个区间上的向量化积分。

希望大家多多支持与鼓励!

论坛优秀回答者

1

主题

9787

帖子

1467

最佳答案
  • 关注者: 334
3#
发表于 2013-10-22 04:17:32 | 只看该作者
总结的不错
另外可以考虑用odesolver求解
而且再比较效率的同时再比较下精度就更加完善了

期待你的下期分享帖

论坛优秀回答者

24

主题

1万

帖子

1638

最佳答案
  • 关注者: 943
4#
 楼主| 发表于 2013-10-22 04:26:10 | 只看该作者

多谢你的建议!上面都是按默认精度在比较,我争取下次分享的时候把精度考虑进去。
另外,能否简要指导一下,odesolver求解方法或思想?

论坛优秀回答者

1

主题

9787

帖子

1467

最佳答案
  • 关注者: 334
5#
发表于 2013-10-22 04:33:40 | 只看该作者
本帖最后由 kaaaf123 于 2013-10-22 04:39 编辑
winner245 发表于 2013-10-22 04:26
多谢你的建议!上面都是按默认精度在比较,我争取下次分享的时候把精度考虑进去。
另外,能否简要指导一 ...

我是说不同积分算法的精度,例如quadv应该是是基于quad,quad的特点是速度快精度低;另外odesolver的精度估计也不是很高

就用你给的例子好了:
t = 10*rand(2,6);
f = @(t)@(x,~) (x.^2-1)./(x.^2+1).*cos(t(:).*x.*exp(-x));
[~,y] = ode45(f(t),[0 .5 1],zeros(size(t)));
y = reshape(y(3,:),size(t));

也可以不用嵌套匿名函数,而直接用参数,中间两句改为
f = @(x,~) (x.^2-1)./(x.^2+1).*cos(t(:).*x.*exp(-x));
[~,y] = ode45(f,[0 .5 1],zeros(size(t)));

新手

22 麦片

财富积分


050


1

主题

69

帖子

1

最佳答案
6#
发表于 2013-10-24 01:05:15 | 只看该作者
This is great.  期待。

新手

30 麦片

财富积分


050


4

主题

56

帖子

2

最佳答案
7#
发表于 2013-10-25 23:50:47 | 只看该作者
头一次听说唉,

新手

10 麦片

财富积分


050


8

主题

72

帖子

0

最佳答案
QQ
8#
发表于 2013-10-26 04:23:17 | 只看该作者
z收藏收藏收藏

论坛优秀回答者

24

主题

1万

帖子

1638

最佳答案
  • 关注者: 943
9#
 楼主| 发表于 2013-10-28 02:35:04 | 只看该作者
本帖最后由 winner245 于 2013-10-28 02:51 编辑
kaaaf123 发表于 2013-10-22 04:33
我是说不同积分算法的精度,例如quadv应该是是基于quad,quad的特点是速度快精度低;另外odesolver的精度 ...

kaaaf123 提供的这个 用 odesolver 求解向量化定积分 的例子实在是太经典了!创意很好、很新颖!不过,kaaaf 兄的代码实在是太精炼了,很多人未必能真正体会到这个方法的精髓。下面,基于我个人的理解加以解释。

1. 先说说 odesolver 求解定积分,尤其是向量化定积分的原理

我们常见的 odesolver 函数,比如:ode45、ode23、ode15s, 都是把方程组里多个微分方程看做一个整体(列向量),通过一次 odesolver 函数调用来求解每个微分方程。如果我们将每个参数值 t 对应的定积分的被积函数看做一个微分方程,我们就得到了一个微分方程组,方程的个数正好等于 t 的不同取值个数。然后,直接用 odesolver 求解该微分方程组,得到的就是我们要求的不同参数 t 对应的积分。

2.  odesolver 函数要求微分方程组是一个函数句柄,且该函数句柄返回的是一个列向量。这就是上述代码中为什么出现了 t(:)

3. odesolver 求解微分方程组时,定义的函数默认是两个自变量,常见的形如: y' = f(x, y),这里的微分是针对第一个变量 x 进行。odesolver 考虑的是一般情形,即微分方程组右边还会出现 y。而本帖考虑的积分全部是 y = int f(x) dx 的情形,也就是微分方程里不包含第二个变量 y。但是,为了符合 odesolver 的语法规则,我们依然需要将 f 定义为两个自变量,第二个自变量可以随意给个名字,也可以像 kaaaf123 那样,直接给个占位符 ~。

4. 正是因为我们的积分函数是 y = f(x) ,而不是 odesolver 的一般情形 y ' = f(x, y),我们可以任意指定 y 在积分起始点 0 处的初值(这是因为 y 并没有出现在 f 中,任意给定的初值并不会改变微分方程组本身)。待我们按任意初值求得积分后,我们还需要用求得的积分数值减去我们指定的初值,这样求得的差值才是我们真正需要的定积分。举个例子,我们将上面代码里的初值改成随机数:
  1. t = 10*rand(2,6);
  2. f = @(t)@(x,~) (x.^2-1)./(x.^2+1).*cos(t(:).*x.*exp(-x));
  3. [~,y] = ode45(f(t),[0 .5 1],rand(size(t)));
  4. y = reshape(y(3,:)-y(1,:),size(t))
复制代码
可以发现,这段代码得到了相同的结果。显然,为了省去最后这一步求差值,最方便的办法莫过于把初值全部设置为0,这也就是 kaaaf123 的做法 (现在,大家能知道上面代码的精妙所在了吧)

5. odesolver 默认是求解积分区间(这里是[0,1])上许多个点的积分,并将每个点的结果全部存储到返回值里,这样,返回值的每一行对应一个点的积分,最后一行就是积分终点1的积分,第一行就是积分初值。但我们这里显然只需要得到最后一个点积分值,中间的积分值是不需要的,为了避免无谓的存储开销,可以将积分区间参数指定为3个变量,即在0~1之间增加一个点(如0.5),变成了 [0 0.5 1],odesolver 就默认只存储这3个点的积分数值了,这样能节省一点内存(但odesolver 内部还是会自动计算每个点的积分,只是不会全部存储,因此,这样做不会对运算效率有明显的影响)。另外,这个点是可以在积分区间里随意选取的,kaaaf123 选的是积分区间的中点0.5。因为我们只关心终点1处的积分,故需要取最后一行(第3行),这就是上面代码里为什么是 y(3,:)

新手

33 麦片

财富积分


050


4

主题

104

帖子

6

最佳答案
10#
发表于 2013-10-28 12:08:47 | 只看该作者
收藏啦!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 快乐扑克3如何跟对子 返回列表