关于一个积分的求教

eryun05 · 发表于 2006-3-20 11:28

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

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

x

有一个二重积分，且内含一个广义积分，单独的二重积分和广义积分程序都可用，但放在一起调用时出错，请各位大侠帮忙。程序在上传的文件里，用fortran编的。被积函数已给出。 
c program d3r9
 common xmax
 pi=3.1415926
 nval=10
 write(*,*)'integral of r^2 over a spherical volume'
 write(*,'(/4x,a,t14,a,t24,a)')'radius','quad3d','actual'
 xmax=nval
 xmin=-xmax
 call quad2d(xmin,xmax,s)
 write(*,'(1x,f8.2,<st1:chmetcnv TCSC="0" NumberType="1" Negative="False" HasSpace="False" SourceValue="2" UnitName="F" w:st="on">2f</st1:chmetcnv>10.4)')xmax,s
 pause
 stop
 end
 
 function func(x,y)
 call qgys(x,y,sh)
 func=(x**2+y**2)*sh
 return
 end
 
 function y1(x,sh)
 common xmax
 y1=-sqrt(abs((xmax**2-x**2)*sh))
 return
 end
 
 function y2(x,sh)
 common xmax
 y2=sqrt(abs((xmax**2-x**2)*sh))
 return
 end
 
 subroutine qgausx(func,a,b,ss)
 dimension x(5),w(5)
 data x/.1488743389,.4333953941,
 * .6794095682,.8650633666,.9739065285/
 data w/.2955242247,.2692667193,
 * .2190863625,.1494513491,.0666713443/
 xm=0.5*(b+a)
 xr=0.5*(b-a)
 ss=0.
 do 10 j=1,5
 dx=xr*x(j)
 ss=ss+w(j)*(func(xm+dx)+func(xm-dx))
10 continue
 ss=xr*ss
 return
 end
 
 subroutine qgausy(func,a,b,ss)
 dimension x(5),w(5)
 data x/.1488743389,.4333953941,
 * .6794095682,.8650633666,.9739065285/
 data w/.2955242247,.2692667193,
 * .2190863625,.1494513491,.0666713443/
 xm=0.5*(b+a)
 xr=0.5*(b-a)
 ss=0.
 do 10 j=1,5
 dx=xr*x(j)
 ss=ss+w(j)*(func(xm+dx)+func(xm-dx))
10 continue
 ss=xr*ss
 return
 end
 
 
 subroutine quad2d(x1,x2,ss)
 external h
 call qgausx(h,x1,x2,ss)
 return
 end
 
 
 function g(yy)
 common /xyz/x,y
 y=yy
 g=func(x,y)
 return
 end
 
 function h(xx)
 external g
 common /xyz/x,y
 x=xx
 call qgausy(g,y1(x,sh),y2(x,sh),ss)
 h=ss
 return
 end
 
 subroutine qgys(x,y,sh)
 EXTERNAL Fun1,fun2
 DIMENSION RI(8)
 NN=0
30 DO 10 I=NN,NN+1 
 AN=2**I 
 RI=SIMPSN(0.0,AN,Fun1,1E-6,1E-2)
10 CONTINUE
 IF (ABS(RI(NN+1)-RI(NN)).LT.1E-7) GOTO 20
 NN=NN+1
 GOTO 30
20 fs2=ri(nn+1)
 fs1=fun2(x,y)
 sh=fs1-fs2
 return
 END
 
 FUNCTION SIMPSN(A,B,Fun1,EPS,H0)
 H=B-A
 RP=Fun1(A)+Fun1(B)
 N=1
2 zkesi=A-0.5*H
 RC=0.0
 DO 4 I=1,N
 zkesi=zkesi+H
4 RC=Fun1(zkesi)+RC
 R2=(RP+4.0*RC)*H/6.0
 IF (N.EQ.1.OR.ABS(H).GT.H0) GOTO 6
 zkesi=R2-R1
 IF (ABS(R2).GE.1.0) zkesi=zkesi/R2
 IF (ABS(zkesi).LT.EPS) GOTO 8
6 R1=R2
 RP=2.0*RC+RP
 N=N+N
 H=0.5*H
 GOTO 2
8 SIMPSN=R2
 RETURN
 END 
 
 FUNCTION Fun1(x,y,zkesi)
 pi=3.1415926
 fun1=exp(-x*cosh(zkesi))*(sin(pi+y)/(cosh(zkesi)-
 * cos(pi+y))+sin(pi-y)/(cosh(zkesi)
 * -cos(pi-y)))/(2*pi)
 
 RETURN
 END
 
 FUNCTION Fun2(x,y)
 pi=3.1415926
 if (y.le.pi.and.y.ge.(-pi)) then
 fun2=exp(x*cos(y))
 else
 fun2=0
 endif
 RETURN
 END 注：被积表达式 (x**2+y**2)*s x**2+y**2=10 而s是一个广义积分，表达式为s=exp(-x*cosh(zkesi))*(sin(pi+y)/(cosh(zkesi)- cos(pi+y))+sin*pi-y)/(cosh(zkesi)-cos(pi-y)))/(2*pi) 积分区间（0，正无穷） <v:shapetype>
</v:shapetype> 
上面的程序一开始是两部分，分别是单独的二重积分和广义积分，且都可用，但放在一起就不可用，请哪位大侠指都。

风花雪月 · 发表于 2006-4-1 07:08

该函数的广义积分在所有的定义域内都成立吗？ 估计你说的一开始好用只是验证了个别数据

[此贴子已经被作者于2006-4-1 7:09:06编辑过]

账号		自动登录	找回密码
密码			我要加入

[其他相关] 关于一个积分的求教

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

回复：（eryun05）关于一个积分的求教