高斯积分中取的数都是一样的,为什么输出的x是不同的?
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
endsubroutine 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
数值问题,需要给出完整代码,完整输入文件。别人才能帮助你判断。
页:
[1]