hao1982 发表于 2006-11-19 12:33

MATLAB语言的高斯积分点计算程序

下面是计算高斯积分点及其权函数的MATLAB子程序,编写程序中只要调用它就可以了
其中,列如你需要计算[-1,1]区间10个高斯点及其权系数,你直接调用
=grule(10)
10个高斯点即存储在XI数组中,对应权系数存储在WI中



function =grule(n)
% =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.
bp=zeros(n,1); wf=bp; 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

ly830129 发表于 2006-11-22 17:08

有没有foortran版的

hao1982 发表于 2006-11-22 19:14

没,自己看看,翻译成FORTRAN即可

wenyue 发表于 2006-11-23 23:19

谢谢，正在编程用高斯数值积分

zwf100 发表于 2008-3-12 17:03

楼主太牛了:loveliness:

皮卡丘 发表于 2011-6-28 14:36

回复 1 # hao1982 的帖子

楼主还在这个群么 可以帮我解决一个问题么

hao1982 发表于 2011-6-30 21:51

回复 7 # 皮卡丘 的帖子

怎么个帮忙法

皮卡丘 发表于 2011-7-1 09:46

谢谢你首先！我试了syms x theta
A=
fun=inline(A,'x','theta')
dblquad(fun,-1,1,0,2*pi)有错误；
试了syms x theta
A=
aa='A'
fun=inline(aa,'x','theta')
dblquad(fun,-1,1,0,2*pi) 也有错误 想问问大侠应该如何修改啊~

liangcheng 发表于 2011-7-4 17:19

问题很详细啊，回去研究研究
！

twb0624 发表于 2012-4-12 13:08

很好的一个程序，不知道好用不。
页: [1]
查看完整版本: MATLAB语言的高斯积分点计算程序