Fortran Coder

标题: 高斯积分中取的数都是一样的,为什么输出的x是不同的? [打印本页]

作者: ylw0710    时间: 2018-1-12 15:26
标题: 高斯积分中取的数都是一样的,为什么输出的x是不同的?
[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


作者: kyra    时间: 2018-1-13 09:09
数值问题,需要给出完整代码,完整输入文件。别人才能帮助你判断。




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2