声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6405|回复: 22

[共享资源] Guass积分点与权系数计算

  [复制链接]
发表于 2011-8-13 15:43 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

Gauss Intergtation.rar

2.47 KB, 下载次数: 38

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2011-8-15 23:59 | 显示全部楼层
这个可以配合
高斯积分点计算程序 http://forum.vibunion.com/thread-32381-1-1.html
是否有人能提供一下高斯积分表 http://forum.vibunion.com/thread-32144-1-1.html
发表于 2011-8-18 13:28 | 显示全部楼层
谢谢楼主
发表于 2011-12-24 23:46 | 显示全部楼层
垃圾~ 浪费了我很多财富,根本下载不下来
发表于 2011-12-24 23:54 | 显示全部楼层
发表于 2011-12-26 23:37 | 显示全部楼层
ChaChing 发表于 2011-12-24 23:54
别生气, 体能还是很容易取得的!

谢谢您.
只是觉得不厚道,下载不了的还在这挂着。
发表于 2011-12-27 00:17 | 显示全部楼层
本帖最后由 bainhome 于 2011-12-27 01:46 编辑

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

Gauss Intergtation.rar

2.47 KB, 下载次数: 2

发表于 2011-12-27 18:04 | 显示全部楼层
bainhome 发表于 2011-12-27 00:17
先看看自己的浏览器下载是不是被迅雷,快车之类给自动嗅探了,很多论坛不支持此类多线程下载,为此专门试了 ...

你不用在这装清高,有可能是迅雷劫持,这个是有可能,但是右键是下载不下来。
发表于 2011-12-27 20:23 | 显示全部楼层
本帖最后由 bainhome 于 2011-12-27 20:23 编辑

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

for sb.jpg
发表于 2011-12-27 22:29 | 显示全部楼层
本帖最后由 ChaChing 于 2011-12-27 22:32 编辑

回复 8 # billy211 的帖子

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

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

点评

很遗憾一个提醒又让老兄回复变得如此尴尬.只是依然无法习惯某些技术上彻底无知、出了问题不假思索责怪他人、同时不知感谢他人努力的喷子.  发表于 2011-12-27 23:07
发表于 2011-12-28 00:43 | 显示全部楼层
回复 9 # bainhome 的帖子

喔, 不觉得尴尬!
基本上, 现在年轻人许多都是如此, 可能个人有点看多了
还有这位同学应该奖励下, 原因:逼高手出来! 你说对否!:@)
发表于 2012-1-3 21:28 | 显示全部楼层
bainhome 发表于 2011-12-27 20:23
首先,我和放附件的楼主以及你都并不认识,所以客观地提醒一下:
我这里都下载并上传了,我想稍微有点脑子 ...

其实你就是一个名副其实的SB,你不是装清高是什么? 你自认为是刻薄,但实际上你是不知脸耻,不成自己是SB,装清高的二B。你牛逼什么?  就是一个装清高的二百五。还厚颜无耻的在这说大道理。
发表于 2012-1-3 21:34 | 显示全部楼层
ChaChing 发表于 2011-12-27 22:29
回复 8 # billy211 的帖子

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

我确实用右键没下载下来,我的浏览器GreenBrowser确实是没下载下来。
但是我不需要这种自命清高,装逼的人来解答。
牛逼什么牛逼!
发表于 2012-1-4 00:02 | 显示全部楼层
本帖最后由 ChaChing 于 2012-1-4 00:06 编辑

回复 13 # billy211 的帖子

这位仁兄, 可别忘记是你先用词不当的!
怎还火气这麼大, 还是该注意下!
以前个人脾气是会删帖禁言的
 楼主| 发表于 2012-1-4 09:40 | 显示全部楼层
没想到给大家添了这么多麻烦!
Gusssfun        ——    进行高斯积分的函数
Guass_Point     ——    区间[a,b]内进行高斯积分的积分点
grule           ——    区间[-1,1]的高斯积分点
G_I             ——    进行高斯积分点调用格式
  1. function Y = Guassfun(u)

  2. Y = [1 + exp(-u)*sin(4*u);
  3.     sin(sqrt(u))];

  4. % Y=[1 + u.^2 - 3 * u.^3 + 4 * u.^5, 1 + u.^2 - 3 * u.^3 ;
  5. %    1 + u.^2 - 3 * u.^3, u.^2 - 3 * u.^3 + 4 * u.^5;
  6. %    ];

  7. end
复制代码
  1. function [bp,wf] = grule(n)

  2. % -------------------------------------------------------------------------
  3. % This code is obtained by "Google" %
  4. % -------------------------------------------------------------------------

  5. %  [bp,wf]=grule(n)
  6. %  This function computes Gauss base points and weight factors
  7. %  using the algorithm given by Davis and Rabinowitz in
  8. %  'Methods of Numerical Integration', page 365, Academic Press, 1975.

  9. %  Input: n - the number of Guass integraltion points   
  10. %  Output: bp - the value of Guass integration points
  11. %          wf - the value of Guass integration weights

  12. bp = zeros(n,1); wf = zeros(n,1);
  13. iter = 2; m = fix((n+1)/2); e1 = n*(n+1);
  14. mm = 4*m-1; t = (pi/(4*n+2))*(3:4:mm); nn = (1-(1-1/n)/(8*n*n));
  15. xo = nn*cos(t);
  16. for j = 1:iter
  17.    pkm1 = 1; pk = xo;
  18.    for k = 2:n
  19.       t1 = xo.*pk;
  20.       pkp1 = t1-pkm1-(t1-pkm1)/k+t1;
  21.       pkm1 = pk; pk = pkp1;
  22.    end
  23.    den = 1.-xo.*xo; d1 = n*(pkm1-xo.*pk); dpn = d1./den;
  24.    d2pn = (2.*xo.*dpn-e1.*pk)./den;
  25.    d3pn = (4*xo.*d2pn+(2-e1).*dpn)./den;
  26.    d4pn = (6*xo.*d3pn+(6-e1).*d2pn)./den;
  27.    u = pk./dpn; v = d2pn./dpn;
  28.    h = -u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
  29.    p = pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
  30.    dp = dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
  31.    h = h-p./dp; xo = xo+h;
  32. end
  33. bp = -xo-h;
  34. fx = d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
  35. wf = 2*(1-bp.^2)./(fx.*fx);
  36. if (m+m) > n, bp(m)=0; end
  37. if ~((m+m) == n), m=m-1; end
  38. jj = 1:m; n1j = (n+1-jj);
  39. bp(n1j) = -bp(jj); wf(n1j) = wf(jj);

  40. % -------------------------------------------------------------------------
  41. end
复制代码
  1. function varginout = Guass_Point(n,a,b)

  2. % -------------------------------------------------------------------------
  3. %
  4. % Guass_Int.m
  5. %
  6. % Calculate the value of Guass integration points and weights for arbitrary interzone [a,b]
  7. %
  8. % Input :    n - the number of Guass integraltion points
  9. %            a - the low boundary of intergral interzone
  10. %            b - the up boundary of intergral interzone
  11. % Output:   varginout{1} = bp - the value of Guass integration points
  12. %           varginout{2} = wf - the value of Guass integration weights
  13. %
  14. % Written by Yan - 21/3/2011
  15. % Contact: yanyongju@ase.buaa.edu.cn
  16. %
  17. % -------------------------------------------------------------------------

  18. % -------------------------------------------------------------------------
  19. %  The Guass integral function can only used within the interzone [-1,1],
  20. %  If we want to get the integration from integral boundary a to b, as [a,b]
  21. %  we can tranform the [a,b] to [-1,1] with the equation
  22. %                a + b       b - a
  23. %           s = -------  +  ------- * t      
  24. %                  2           2
  25. %  where t is in the interzone [-1,1], while s in [a,b]

  26. %  so we can get
  27. %                   / a + b     b - a     \   b - a
  28. %        f(s)ds = f| ------- + ------- * t | -------dt
  29. %                   \   2         2       /     2
  30. % -------------------------------------------------------------------------

  31. if nargin == 1
  32.     a = -1;
  33.     b = 1;
  34. end

  35. [bp,wf] = grule(n);

  36. sp = size(bp);

  37. gpoint = (a+b)/2*ones(sp) + (b-a)/2 * bp;
  38. wpoint = (b-a)/2 * wf;

  39. varginout = {gpoint, wpoint};

  40. % -------------------------------------------------------------------------
  41. 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. par = Guass_Point(N,a,b);

  21. point = par{1}; weight = par{2};

  22. % -------------------------------------------------------------------------
  23. % func = Guassfun(point);
  24. % % func = feval('Guassfun',point);
  25. %
  26. % [m,n] = size(func);
  27. %
  28. % for i = 1:m
  29. %     for j = 1:n
  30. %         val(i,j) = dot(func{i,j},weight);
  31. %     end
  32. % end
  33. % -------------------------------------------------------------------------
  34. val = 0.0;
  35. for k = 1:N   
  36.     func = Guassfun(point(k));
  37.     val = val + func * weight(k);   
  38. end
  39. % -------------------------------------------------------------------------
  40. end
复制代码

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-5 02:37 , Processed in 0.094905 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表