Guass积分点与权系数计算
【与大家分享一下】在计算过程中常用到Guass积分点的计算。程序grule(n)，用于计算区间[-1,1]上的高斯点及其相应权系数。如果上下限为,积分变量为x，作变换为x=(a+b)/2+t*(b-a)/2，并在此基础上写了适用于任何区间的Guass积分计算程序，详见附件。这个可以配合
高斯积分点计算程序 http://forum.vibunion.com/thread-32381-1-1.html
是否有人能提供一下高斯积分表 http://forum.vibunion.com/thread-32144-1-1.html
谢谢楼主{:{39}:} 垃圾~ 浪费了我很多财富,根本下载不下来 billy211 发表于 2011-12-24 23:46 static/image/common/back.gif
垃圾~ 浪费了我很多财富,根本下载不下来
别生气, 体能还是很容易取得的! ChaChing 发表于 2011-12-24 23:54 static/image/common/back.gif
别生气, 体能还是很容易取得的!
谢谢您.
只是觉得不厚道,下载不了的还在这挂着。 本帖最后由 bainhome 于 2011-12-27 01:46 编辑
先看看自己的浏览器下载是不是被迅雷，快车之类给自动嗅探了，很多论坛不支持此类多线程下载，为此专门试了试：我这里下载正常。如果实在不会解除迅雷之类对页面资源的嗅探，不妨用【右键】该资源→【另存为】试试。
如果万一，我说万一...真是这样令人惊讶的低级错误，“垃圾”、“不厚道”之类说法可就有点丢人了哦。
bainhome 发表于 2011-12-27 00:17 static/image/common/back.gif
先看看自己的浏览器下载是不是被迅雷，快车之类给自动嗅探了，很多论坛不支持此类多线程下载，为此专门试了 ...
你不用在这装清高,有可能是迅雷劫持,这个是有可能,但是右键是下载不下来。
本帖最后由 bainhome 于 2011-12-27 20:23 编辑
首先，我和放附件的楼主以及你都并不认识，所以客观地提醒一下：
我这里都下载并上传了，我想稍微有点脑子的都会意识到是自己电脑的问题。居然这么明显的情况，还死不要脸嘴硬，S*B到这个程度，只好给你无限真诚的祝福了。
本帖最后由 ChaChing 于 2011-12-27 22:32 编辑
回复 8 # billy211 的帖子
汗, 我是不清楚"被迅雷，快车"是什麼意思? 又懒得google:@L
不过我使用ie点选或右键下载是没问题的! 其实其中主程序grule在2F的连接即有
最后别人的好意帮忙, 没感谢便罢, 也别当成...:@)
回复 9 # bainhome 的帖子
喔, 不觉得尴尬!
基本上, 现在年轻人许多都是如此, 可能个人有点看多了
还有这位同学应该奖励下, 原因:逼高手出来! 你说对否!:@) bainhome 发表于 2011-12-27 20:23 static/image/common/back.gif
首先，我和放附件的楼主以及你都并不认识，所以客观地提醒一下：
我这里都下载并上传了，我想稍微有点脑子 ...
其实你就是一个名副其实的SB，你不是装清高是什么? 你自认为是刻薄，但实际上你是不知脸耻，不成自己是SB，装清高的二B。你牛逼什么?就是一个装清高的二百五。还厚颜无耻的在这说大道理。 ChaChing 发表于 2011-12-27 22:29 static/image/common/back.gif
回复 8 # billy211 的帖子
汗, 我是不清楚"被迅雷，快车"是什麼意思? 又懒得google
我确实用右键没下载下来，我的浏览器GreenBrowser确实是没下载下来。
但是我不需要这种自命清高，装逼的人来解答。
牛逼什么牛逼！ 本帖最后由 ChaChing 于 2012-1-4 00:06 编辑
回复 13 # billy211 的帖子
这位仁兄, 可别忘记是你先用词不当的!
怎还火气这麼大, 还是该注意下!
以前个人脾气是会删帖禁言的
没想到给大家添了这么多麻烦！
Gusssfun —— 进行高斯积分的函数
Guass_Point —— 区间内进行高斯积分的积分点
grule —— 区间[-1,1]的高斯积分点
G_I —— 进行高斯积分点调用格式function Y = Guassfun(u)
Y = [1 + exp(-u)*sin(4*u);
sin(sqrt(u))];
% Y=[1 + u.^2 - 3 * u.^3 + 4 * u.^5, 1 + u.^2 - 3 * u.^3 ;
% 1 + u.^2 - 3 * u.^3, u.^2 - 3 * u.^3 + 4 * u.^5;
% ];
endfunction = grule(n)
% -------------------------------------------------------------------------
% This code is obtained by "Google" %
% -------------------------------------------------------------------------
%=grule(n)
%This function computes Gauss base points and weight factors
%using the algorithm given by Davis and Rabinowitz in
%'Methods of Numerical Integration', page 365, Academic Press, 1975.
%Input: n - the number of Guass integraltion points
%Output: bp - the value of Guass integration points
% wf - the value of Guass integration weights
bp = zeros(n,1); wf = zeros(n,1);
iter = 2; m = fix((n+1)/2); e1 = n*(n+1);
mm = 4*m-1; t = (pi/(4*n+2))*(3:4:mm); nn = (1-(1-1/n)/(8*n*n));
xo = nn*cos(t);
for j = 1:iter
pkm1 = 1; pk = xo;
for k = 2:n
t1 = xo.*pk;
pkp1 = t1-pkm1-(t1-pkm1)/k+t1;
pkm1 = pk; pk = pkp1;
end
den = 1.-xo.*xo; d1 = n*(pkm1-xo.*pk); dpn = d1./den;
d2pn = (2.*xo.*dpn-e1.*pk)./den;
d3pn = (4*xo.*d2pn+(2-e1).*dpn)./den;
d4pn = (6*xo.*d3pn+(6-e1).*d2pn)./den;
u = pk./dpn; v = d2pn./dpn;
h = -u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
p = pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
dp = dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
h = h-p./dp; xo = xo+h;
end
bp = -xo-h;
fx = d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
wf = 2*(1-bp.^2)./(fx.*fx);
if (m+m) > n, bp(m)=0; end
if ~((m+m) == n), m=m-1; end
jj = 1:m; n1j = (n+1-jj);
bp(n1j) = -bp(jj); wf(n1j) = wf(jj);
% -------------------------------------------------------------------------
end
function varginout = Guass_Point(n,a,b)
% -------------------------------------------------------------------------
%
% Guass_Int.m
%
% Calculate the value of Guass integration points and weights for arbitrary interzone
%
% Input : n - the number of Guass integraltion points
% a - the low boundary of intergral interzone
% b - the up boundary of intergral interzone
% Output: varginout{1} = bp - the value of Guass integration points
% varginout{2} = wf - the value of Guass integration weights
%
% Written by Yan - 21/3/2011
% Contact: yanyongju@ase.buaa.edu.cn
%
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%The Guass integral function can only used within the interzone [-1,1],
%If we want to get the integration from integral boundary a to b, as
%we can tranform the to [-1,1] with the equation
% a + b b - a
% s = -------+------- * t
% 2 2
%where t is in the interzone [-1,1], while s in
%so we can get
% / a + b b - a \ b - a
% f(s)ds = f| ------- + ------- * t | -------dt
% \ 2 2 / 2
% -------------------------------------------------------------------------
if nargin == 1
a = -1;
b = 1;
end
= grule(n);
sp = size(bp);
gpoint = (a+b)/2*ones(sp) + (b-a)/2 * bp;
wpoint = (b-a)/2 * wf;
varginout = {gpoint, wpoint};
% -------------------------------------------------------------------------
end
function val = G_I(a,b,N)
% a = 0; b = 1; N = 10;
%
% -------------------------------------------------------------------------
%
% Guass_Int.m
%
% Calculate the value of @Guassfun by Guass intergraltion
%
% Input :fun - the function of Guass integraltion
% N - the number of Guass integraltion points
% a - the low boundary of intergral interzone
% b - the up boundary of intergral interzone
% Output:val - the value of Guass integraltion
%
% Written by Yan - 21/3/2011
% Contact: yanyongju@ase.buaa.edu.cn
%
% -------------------------------------------------------------------------
par = Guass_Point(N,a,b);
point = par{1}; weight = par{2};
% -------------------------------------------------------------------------
% func = Guassfun(point);
% % func = feval('Guassfun',point);
%
% = size(func);
%
% for i = 1:m
% for j = 1:n
% val(i,j) = dot(func{i,j},weight);
% end
% end
% -------------------------------------------------------------------------
val = 0.0;
for k = 1:N
func = Guassfun(point(k));
val = val + func * weight(k);
end
% -------------------------------------------------------------------------
end
页:
[1]
2