Fortran Coder

查看: 202|回复: 1

[微积分] 高斯积分中取的数都是一样的,为什么输出的x是不同的?

[复制链接]

12

帖子

8

主题

0

精华

入门

F 币
89 元
贡献
55 点
发表于 2018-1-12 15:26:32 | 显示全部楼层 |阅读模式
[Fortran] 纯文本查看 复制代码
subroutine LGAUS220(a,b,f,eps,g,Xf,xm,AAA,zzz,n1,nucl)
	external f0
      dimension t(5),c(5)
      double precision a,b,f11,g,t,c,s,p,h,aa,w,x,q,X1,AAA,eps,xw,f12,
     *	xn,xm,f,zzz,bb,f0,xf
      integer n,kk,iset,n1,nucl
	data t/-0.9061798459,-0.5384693101,0.0,
     &       0.5384693101,0.9061798459/
      data c/0.2369268851,0.4786286705,0.568888889,
     &       0.4786286705,0.2369268851/
C      DATA t/ -0.993128599185094924786,       -0.963971927277913791268,
C     &	    -0.912234428251325905868,       -0.839116971822218823395,
C     &	    -0.746331906460150192617,       -0.636053680726515025453,
C     &	    -0.510867001950827098004,       -0.373706088715419560673,
C     &	    -0.227785851141645078080,       -0.076526521133497333755,
C     &         0.076526521133497333755,        0.227785851141645078080,
C     &	     0.373706088715419560673,        0.510867001950827098004,
C     &         0.636053680726515025453,        0.746331906460150192617,
C     &	     0.839116971822218823395,        0.912234428251325905868,
C     &         0.963971927277913791268,        0.993128599185094924786/
C	DATA c/  0.017614007139152118312,        0.040601429800386941331,
C     &	     0.062672048334109063570,        0.083276741576704748725,
C     &	     0.101930119817240435037,        0.118194531961518417312,
C     &         0.131688638449176626898,        0.142096109318382051329,
C     &         0.149172986472603746788,        0.152753387130725850698,
C     &         0.152753387130725850698,        0.149172986472603746788,
C     &         0.142096109318382051329,        0.131688638449176626898,
C     &         0.118194531961518417312,        0.101930119817240435037,
C     &         0.083276741576704748725,        0.062672048334109063570,
C     &         0.040601429800386941331,        0.017614007139152118312/
       
            
	
      m=1
      s=(b-a)*0.001
      p=0.0
10    h=(b-a)/m
      g=0.0
      do 30 i=1,m
         aa=a+(i-1)*h
         bb=a+i*h
         w=0.0
	   do 20 kk=1,5
         x=((bb-aa)*t(kk)+(bb+aa))/2.0
        
         w=w+f0(n1,x,xf,xM,Aaa,zzz,nucl)*c(kk)
      
20       continue
         g=g+w

30    continue
      g=g*h/2.0
	
      q=abs(g-p)/(1.0+abs(g))
      
      if((q.ge.eps).and.(abs(h).gt.abs(s)))then
        p=g
        m=m+1
      
        goto 10
	end if
	return
      end
[Fortran] 纯文本查看 复制代码
subroutine LGAUS20(a,b,f,eps,g,Xf,xm,AAA,zzz,nucl,n1,bb0)
	external f1
      dimension t(5),c(5)
      double precision a,b,f11,g,t,c,s,p,h,aa,bb,w,x,q,X1,AAA,eps,xw,f12
     *	,xn,xm,f,zzz,f1,bb0,xf
      integer n,kk,iset,n1,nucl
	data t/-0.9061798459,-0.5384693101,0.0,
     &       0.5384693101,0.9061798459/
      data c/0.2369268851,0.4786286705,0.568888889,
     &       0.4786286705,0.2369268851/
C      DATA t/ -0.993128599185094924786,       -0.963971927277913791268,
C     &	    -0.912234428251325905868,       -0.839116971822218823395,
C     &	    -0.746331906460150192617,       -0.636053680726515025453,
C     &	    -0.510867001950827098004,       -0.373706088715419560673,
C     &	    -0.227785851141645078080,       -0.076526521133497333755,
C     &         0.076526521133497333755,        0.227785851141645078080,
C     &	     0.373706088715419560673,        0.510867001950827098004,
C     &         0.636053680726515025453,        0.746331906460150192617,
C     &	     0.839116971822218823395,        0.912234428251325905868,
C     &         0.963971927277913791268,        0.993128599185094924786/
C	DATA c/  0.017614007139152118312,        0.040601429800386941331,
C     &	     0.062672048334109063570,        0.083276741576704748725,
C     &	     0.101930119817240435037,        0.118194531961518417312,
C     &         0.131688638449176626898,        0.142096109318382051329,
C     &         0.149172986472603746788,        0.152753387130725850698,
C     &         0.152753387130725850698,        0.149172986472603746788,
C     &         0.142096109318382051329,        0.131688638449176626898,
C     &         0.118194531961518417312,        0.101930119817240435037,
C     &         0.083276741576704748725,        0.062672048334109063570,
C     &         0.040601429800386941331,        0.017614007139152118312/
       
       
      
      
	
      m=1
      s=(b-a)*0.001
      p=0.0
10    h=(b-a)/m
      g=0.0
      do 30 i=1,m
         aa=a+(i-1)*h
         bb=a+i*h
         w=0.0
	   do 20 kk=1,5
         x=((bb-aa)*t(kk)+(bb+aa))/2.0
      
         w=w+f1(n1,bb0,x,xm,Xf,Aaa,Zzz,nucl)*c(kk)
      
20       continue
	
         g=g+w

30    continue
      g=g*h/2.0
	
      q=abs(g-p)/(1.0+abs(g))
      
      if((q.ge.eps).and.(abs(h).gt.abs(s)))then
        p=g
        m=m+1
      
        goto 10
	end if
	return
      end

回复

使用道具 举报

71

帖子

0

主题

0

精华

熟手

F 币
228 元
贡献
106 点

新人勋章美女勋章元老勋章热心勋章规矩勋章

QQ
发表于 2018-1-13 09:09:50 | 显示全部楼层
数值问题,需要给出完整代码,完整输入文件。别人才能帮助你判断。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2018-4-22 18:24

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

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