求解：高斯－雷让德积分点

rita

用matlab的什么命令可以求解n个高斯－雷让德积分点

rita

用什么命令可以得出高斯－勒让德求积公式的节点

w89986581

function [bp,wf]=gausspoint(n)
%%%计算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);

rita

怎么用呀。我现在是划分节点，采用高斯－勒让德积分点作为节点，请问怎么办到

w89986581

[bp,wf]=gausspoint(n)中的n表示阶数，bp为节点。若你的上下限不是[-1,1]，就把它映射到[-1,1]即可。
参考李庆扬, 王能超, 易大义. 数值分析[M]. 湖北: 华中科技大学出版社, 1986

rita

matlab找不到这个gausspoint函数呀

w89986581

:@L 这不是个现成的gausspoint程序么？
把它保存到m起来，文件名和函数一样以后，就可以调用了.....无语了

rita

知道了，谢谢

[ 本帖最后由 eight 于 2007-6-8 15:21 编辑 ]

mxlzhenzhu

w89986581 发表于 2007-6-8 14:23

                               
登录/注册后可看大图

=gausspoint(n)中的n表示阶数，bp为节点。若你的上下限不是[-1,1]，就把它映射到[-1,1]即可。
参考李庆扬, ...

您好！我的积分对象是一个三维六面体单元，那么怎么确定积分上下限的函数表达式呢？ 需要最小二乘拟合平面出来么？我有点儿疑惑。