查看: 31872|回復: 56|關注: 0
打印 上一主題 下一主題

山东体彩快乐扑克三: MATLAB高效編程之向量化積分

  [復制鏈接]

論壇優秀回答者

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

主題

1萬

帖子

1639

最佳答案
  • 關注者: 972
跳轉到指定樓層
1#
發表于 2013-10-22 03:30:23 | 只看該作者 回帖獎勵 |倒序瀏覽 |閱讀模式
本帖最后由 winner245 于 2013-10-24 01:54 編輯

前些天和吳兄(編著《MATLAB高效編程技巧與應用:25個案例分析》)閑談中,談到了如何提高數值積分的效率,達成的共識是利用向量化積分來實現高效率的數值積分。

眾所周知,MATLAB為我們提供了兩個非常易用的向量化積分函數:quadv 和 integral 函數。本帖的目的是探討一下,如何利用這兩個積分函數,解決三類不同的向量化積分問題。

我先拋磚引玉,在跟帖里給出我個人這些天的體會?;隊蠹也斡胩致?,批評指正,希望能把這個話題繼續下去?;隊魏問禱值幕疤?/font>

論壇優秀回答者

24

主題

1萬

帖子

1639

最佳答案
  • 關注者: 972
來自 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

帖子

1470

最佳答案
  • 關注者: 341
3#
發表于 2013-10-22 04:17:32 | 只看該作者
總結的不錯
另外可以考慮用odesolver求解
而且再比較效率的同時再比較下精度就更加完善了

期待你的下期分享帖

論壇優秀回答者

24

主題

1萬

帖子

1639

最佳答案
  • 關注者: 972
4#
 樓主| 發表于 2013-10-22 04:26:10 | 只看該作者

多謝你的建議!上面都是按默認精度在比較,我爭取下次分享的時候把精度考慮進去。
另外,能否簡要指導一下,odesolver求解方法或思想?

論壇優秀回答者

1

主題

9787

帖子

1470

最佳答案
  • 關注者: 341
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)));

新手

23 麥片

財富積分


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萬

帖子

1639

最佳答案
  • 關注者: 972
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 | 只看該作者
收藏啦!
您需要登錄后才可以回帖 登錄 | 注冊

本版積分規則

關閉

站長推薦上一條 /4 下一條

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