声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6707|回复: 9

[编程技巧] MATLAB语言的高斯积分点计算程序

[复制链接]
发表于 2006-11-19 12:33 | 显示全部楼层 |阅读模式

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

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

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



function [bp,wf]=grule(n)
% [bp,wf]=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

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2006-11-22 17:08 | 显示全部楼层

有没有foortran版的

 楼主| 发表于 2006-11-22 19:14 | 显示全部楼层

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

发表于 2006-11-23 23:19 | 显示全部楼层
谢谢,正在编程用高斯数值积分
发表于 2008-3-12 17:03 | 显示全部楼层
楼主太牛了:loveliness:
发表于 2011-6-28 14:36 | 显示全部楼层
回复 1 # hao1982 的帖子

楼主还在这个群么 可以帮我解决一个问题么
 楼主| 发表于 2011-6-30 21:51 | 显示全部楼层
回复 7 # 皮卡丘 的帖子

怎么个帮忙法
发表于 2011-7-1 09:46 | 显示全部楼层
谢谢你首先!我试了syms x theta
A=[x*theta,cos(theta);x*sin(theta),x^2+cos(theta)]
fun=inline(A,'x','theta')
dblquad(fun,-1,1,0,2*pi)有错误;
试了syms x theta
A=[x*theta,cos(theta);x*sin(theta),x^2+cos(theta)]
aa='A'
fun=inline(aa,'x','theta')
dblquad(fun,-1,1,0,2*pi) 也有错误 想问问大侠应该如何修改啊~
发表于 2011-7-4 17:19 | 显示全部楼层
问题很详细啊,回去研究研究
发表于 2012-4-12 13:08 | 显示全部楼层
很好的一个程序,不知道好用不。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-5 11:39 , Processed in 0.068650 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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