Guass积分点与权系数计算

yanyongju

【与大家分享一下】在计算过程中常用到Guass积分点的计算。程序grule(n)，用于计算区间[-1,1]上的高斯点及其相应权系数。如果上下限为[a,b],积分变量为x，作变换为x=(a+b)/2+t*(b-a)/2，并在此基础上写了适用于任何区间的Guass积分计算程序，详见附件。

ChaChing

这个可以配合
高斯积分点计算程序 http://forum.vibunion.com/thread-32381-1-1.html
是否有人能提供一下高斯积分表 http://forum.vibunion.com/thread-32144-1-1.html

lordgo

谢谢楼主

billy211

垃圾~ 浪费了我很多财富,根本下载不下来

ChaChing

billy211 发表于 2011-12-24 23:46

                               
登录/注册后可看大图

垃圾~ 浪费了我很多财富,根本下载不下来

别生气, 体能还是很容易取得的!

billy211

ChaChing 发表于 2011-12-24 23:54

                               
登录/注册后可看大图

别生气, 体能还是很容易取得的!

谢谢您.
只是觉得不厚道,下载不了的还在这挂着。

bainhome

先看看自己的浏览器下载是不是被迅雷，快车之类给自动嗅探了，很多论坛不支持此类多线程下载，为此专门试了试：我这里下载正常。如果实在不会解除迅雷之类对页面资源的嗅探，不妨用【右键】该资源→【另存为】试试。
如果万一，我说万一...真是这样令人惊讶的低级错误，“垃圾”、“不厚道”之类说法可就有点丢人了哦。

billy211

bainhome 发表于 2011-12-27 00:17

                               
登录/注册后可看大图

先看看自己的浏览器下载是不是被迅雷，快车之类给自动嗅探了，很多论坛不支持此类多线程下载，为此专门试了 ...

你不用在这装清高,有可能是迅雷劫持,这个是有可能,但是右键是下载不下来。

bainhome

首先，我和放附件的楼主以及你都并不认识，所以客观地提醒一下：
我这里都下载并上传了，我想稍微有点脑子的都会意识到是自己电脑的问题。居然这么明显的情况，还死不要脸嘴硬，S*B到这个程度，只好给你无限真诚的祝福了。

ChaChing

回复 8 # billy211 的帖子

汗, 我是不清楚"被迅雷，快车"是什麼意思? 又懒得google:@L
不过我使用ie点选或右键下载是没问题的! 其实其中主程序grule在2F的连接即有

最后别人的好意帮忙, 没感谢便罢, 也别当成...:@)

ChaChing

回复 9 # bainhome 的帖子

喔, 不觉得尴尬!
基本上, 现在年轻人许多都是如此, 可能个人有点看多了
还有这位同学应该奖励下, 原因:逼高手出来! 你说对否!:@)

billy211

bainhome 发表于 2011-12-27 20:23

                               
登录/注册后可看大图

首先，我和放附件的楼主以及你都并不认识，所以客观地提醒一下：
我这里都下载并上传了，我想稍微有点脑子 ...

其实你就是一个名副其实的SB，你不是装清高是什么? 你自认为是刻薄，但实际上你是不知脸耻，不成自己是SB，装清高的二B。你牛逼什么?  就是一个装清高的二百五。还厚颜无耻的在这说大道理。

billy211

ChaChing 发表于 2011-12-27 22:29

                               
登录/注册后可看大图

回复 8 # billy211 的帖子

汗, 我是不清楚"被迅雷，快车"是什麼意思? 又懒得google

我确实用右键没下载下来，我的浏览器GreenBrowser确实是没下载下来。
但是我不需要这种自命清高，装逼的人来解答。
牛逼什么牛逼！

ChaChing

回复 13 # billy211 的帖子

这位仁兄, 可别忘记是你先用词不当的!
怎还火气这麼大, 还是该注意下!
以前个人脾气是会删帖禁言的

yanyongju

没想到给大家添了这么多麻烦！
Gusssfun        ——    进行高斯积分的函数
Guass_Point     ——    区间[a,b]内进行高斯积分的积分点
grule           ——    区间[-1,1]的高斯积分点
G_I             ——    进行高斯积分点调用格式

1. function Y = Guassfun(u)
   
2. 
   
3. Y = [1 + exp(-u)*sin(4*u);
   
4.     sin(sqrt(u))];
   
5. 
   
6. % Y=[1 + u.^2 - 3 * u.^3 + 4 * u.^5, 1 + u.^2 - 3 * u.^3 ;
   
7. %    1 + u.^2 - 3 * u.^3, u.^2 - 3 * u.^3 + 4 * u.^5;
   
8. %    ];
   
9. 
   
10. end

复制代码

1. function [bp,wf] = grule(n)
   
2. 
   
3. % -------------------------------------------------------------------------
   
4. % This code is obtained by "Google" %
   
5. % -------------------------------------------------------------------------
   
6. 
   
7. %  [bp,wf]=grule(n)
   
8. %  This function computes Gauss base points and weight factors
   
9. %  using the algorithm given by Davis and Rabinowitz in
   
10. %  'Methods of Numerical Integration', page 365, Academic Press, 1975.
    
11. 
    
12. %  Input: n - the number of Guass integraltion points   
    
13. %  Output: bp - the value of Guass integration points
    
14. %          wf - the value of Guass integration weights
    
15. 
    
16. bp = zeros(n,1); wf = zeros(n,1);
    
17. iter = 2; m = fix((n+1)/2); e1 = n*(n+1);
    
18. mm = 4*m-1; t = (pi/(4*n+2))*(3:4:mm); nn = (1-(1-1/n)/(8*n*n));
    
19. xo = nn*cos(t);
    
20. for j = 1:iter
    
21.    pkm1 = 1; pk = xo;
    
22.    for k = 2:n
    
23.       t1 = xo.*pk;
    
24.       pkp1 = t1-pkm1-(t1-pkm1)/k+t1;
    
25.       pkm1 = pk; pk = pkp1;
    
26.    end
    
27.    den = 1.-xo.*xo; d1 = n*(pkm1-xo.*pk); dpn = d1./den;
    
28.    d2pn = (2.*xo.*dpn-e1.*pk)./den;
    
29.    d3pn = (4*xo.*d2pn+(2-e1).*dpn)./den;
    
30.    d4pn = (6*xo.*d3pn+(6-e1).*d2pn)./den;
    
31.    u = pk./dpn; v = d2pn./dpn;
    
32.    h = -u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
    
33.    p = pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
    
34.    dp = dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
    
35.    h = h-p./dp; xo = xo+h;
    
36. end
    
37. bp = -xo-h;
    
38. fx = d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
    
39. wf = 2*(1-bp.^2)./(fx.*fx);
    
40. if (m+m) > n, bp(m)=0; end
    
41. if ~((m+m) == n), m=m-1; end
    
42. jj = 1:m; n1j = (n+1-jj);
    
43. bp(n1j) = -bp(jj); wf(n1j) = wf(jj);
    
44. 
    
45. % -------------------------------------------------------------------------
    
46. end
    

复制代码

1. function varginout = Guass_Point(n,a,b)
   
2. 
   
3. % -------------------------------------------------------------------------
   
4. %
   
5. % Guass_Int.m
   
6. %
   
7. % Calculate the value of Guass integration points and weights for arbitrary interzone [a,b]
   
8. %
   
9. % Input :    n - the number of Guass integraltion points
   
10. %            a - the low boundary of intergral interzone
    
11. %            b - the up boundary of intergral interzone
    
12. % Output:   varginout{1} = bp - the value of Guass integration points
    
13. %           varginout{2} = wf - the value of Guass integration weights
    
14. %
    
15. % Written by Yan - 21/3/2011
    
16. % Contact: yanyongju@ase.buaa.edu.cn
    
17. %
    
18. % -------------------------------------------------------------------------
    
19. 
    
20. % -------------------------------------------------------------------------
    
21. %  The Guass integral function can only used within the interzone [-1,1],
    
22. %  If we want to get the integration from integral boundary a to b, as [a,b]
    
23. %  we can tranform the [a,b] to [-1,1] with the equation
    
24. %                a + b       b - a
    
25. %           s = -------  +  ------- * t      
    
26. %                  2           2
    
27. %  where t is in the interzone [-1,1], while s in [a,b]
    
28. 
    
29. %  so we can get
    
30. %                   / a + b     b - a     \   b - a
    
31. %        f(s)ds = f| ------- + ------- * t | -------dt
    
32. %                   \   2         2       /     2
    
33. % -------------------------------------------------------------------------
    
34. 
    
35. if nargin == 1
    
36.     a = -1;
    
37.     b = 1;
    
38. end
    
39. 
    
40. [bp,wf] = grule(n);
    
41. 
    
42. sp = size(bp);
    
43. 
    
44. gpoint = (a+b)/2*ones(sp) + (b-a)/2 * bp;
    
45. wpoint = (b-a)/2 * wf;
    
46. 
    
47. varginout = {gpoint, wpoint};
    
48. 
    
49. % -------------------------------------------------------------------------
    
50. end
    

复制代码

1. function val = G_I(a,b,N)
   
2. % a = 0; b = 1; N = 10;
   
3. %
   
4. % -------------------------------------------------------------------------
   
5. %
   
6. % Guass_Int.m
   
7. %
   
8. % Calculate the value of @Guassfun by Guass intergraltion
   
9. %
   
10. % Input :  fun - the function of Guass integraltion
    
11. %            N - the number of Guass integraltion points
    
12. %            a - the low boundary of intergral interzone
    
13. %            b - the up boundary of intergral interzone
    
14. % Output:  val - the value of Guass integraltion
    
15. %
    
16. % Written by Yan - 21/3/2011
    
17. % Contact: yanyongju@ase.buaa.edu.cn
    
18. %
    
19. % -------------------------------------------------------------------------
    
20. 
    
21. par = Guass_Point(N,a,b);
    
22. 
    
23. point = par{1}; weight = par{2};
    
24. 
    
25. % -------------------------------------------------------------------------
    
26. % func = Guassfun(point);
    
27. % % func = feval('Guassfun',point);
    
28. %
    
29. % [m,n] = size(func);
    
30. %
    
31. % for i = 1:m
    
32. %     for j = 1:n
    
33. %         val(i,j) = dot(func{i,j},weight);
    
34. %     end
    
35. % end
    
36. % -------------------------------------------------------------------------
    
37. val = 0.0;
    
38. for k = 1:N   
    
39.     func = Guassfun(point(k));
    
40.     val = val + func * weight(k);   
    
41. end
    
42. % -------------------------------------------------------------------------
    
43. end
    

复制代码