主题:贝塔PDF中的积分表达式
!========================存在奇异点的积分===================================
subroutine Ymean1(alpha0,beta0,k1,i1,j1)
use MODGLOB
use IMSL
USE INTEGRAL
implicit none
real*8 alpha0,beta0,es,sum0,sum01,sum02,sum03,width1,width2,width3,sum11,sum12,sum13,sum1,width11,width12,width13,&
& sum3,sum00,sum10,width
INTEGER,PARAMETER :: N=81
real*8 f(12),x1(41),x2(20),x3(20),fun1(41),fun2(20),fun3(20),jf1(11),jf11(11)
real*8 f3(12),x31(41),x32(20),x33(20),fun31(41),fun32(20),fun33(20)
real*8,external:: func,fun
real datas(n)
integer ii,jj,k1,i1,j1
es=1.0d-6
f=(/es,1.0d-5,1.0d-4,1.0d-3,1.0d-2,0.1d0,0.9d0,1.0-1.0d-2,1.0-1.0d-3,1.0-1.0d-4,1.0-1.0d-5,1.0-es/)
do ii=1,12
call enlg3(f0,yy(1:81,k1,i1,j1),81,f(ii),f3(ii))
enddo
sum0=es**alpha0/alpha0+es**beta0/beta0
sum1=es**alpha0/alpha0*yy(1,k1,i1,j1)+es**beta0/beta0*yy(81,k1,i1,j1)
sum00=0.0
sum10=0.0
if(ii==6) then
width1=(f(ii+1)-f(ii))/40.d0
width11=(f3(ii+1)-f3(ii))/40.d0
do jj=1,41
x1(jj)=f(ii)+width1*(jj-1)
x31(jj)=f3(ii)+width11*(jj-1)
fun1(jj)=x1(jj)**(alpha0-1.d0)*(1.d0-x1(jj))**(beta0-1.d0)
fun31(jj)=x31(jj)*x1(jj)**(alpha0-1.d0)*(1.d0-x1(jj))**(beta0-1.d0)
enddo
call GenerateData(DATAS,WIDTH,FUNC)
jf1(6)=Trape_Integral(fun1,width1)
jf11(6)=Trape_Integral(fun31,width1)
sum01=sum00+jf1(6)
sum11=sum10+jf11(6)
else if(ii<=5) then
width2=(f(ii+1)-f(ii))/20.d0
width12=(f3(ii+1)-f3(ii))/20.d0
do jj=1,20
x2(jj)=f(ii)+width2*(jj-1)
x32(jj)=f3(ii)+width12*(jj-1)
fun2(jj)=x2(jj)**(alpha0-1.0)*(1.d0-x2(jj))**(beta0-1.d0)
fun32(jj)=x32(jj)*x2(jj)**(alpha0-1.0)*(1.d0-x2(jj))**(beta0-1.d0)
enddo
do ii=1,5
jf1(ii)=Trape_Integral(fun2,width2)
jf11(ii)=Trape_Integral(fun32,width2)
sum02=sum00+jf1(ii)
sum12=sum10+jf11(ii)
enddo
else
width3=(f(ii+1)-f(ii))/20.d0
width13=(f3(ii+1)-f3(ii))/20.d0
do jj=1,20
x3(jj)=f(ii)+width3*(jj-1)
x33(jj)=f3(ii)+width13*(jj-1)
fun3(jj)=x3(jj)**(alpha0-1.0)*(1.d0-x3(jj))**(beta0-1.d0)
fun33(jj)=x33(jj)*x3(jj)**(alpha0-1.0)*(1.d0-x3(jj))**(beta0-1.d0)
enddo
do ii=6,11
jf1(ii)=Trape_Integral(fun3,width3)
jf11(ii)=Trape_Integral(fun33,width3)
sum01=sum00+jf1(ii)
sum11=sum10+jf11(ii)
enddo
endif
sum0=sum0+sum01+sum02+sum03
sum1=sum1+sum11+sum12+sum13
YYM(k1,i1,j1)=sum1/sum0
return
end subroutine
!=======================无奇异点的积分====================================
subroutine Ymean2(k1,i1,j1)
use MODGLOB
use IMSL
USE INTEGRAL
implicit none
real*8 sum,width1,width2,width
INTEGER,PARAMETER :: N=81
real datas(N)
real*8 x(11),fun(11),y1(11),p(11),jf(11)
real*8 ,external:: func
integer ii,jj,k1,i1,j1
sum=0.d0
do ii=1,80
width=(f0(ii+1)-f0(ii))/10.d0
width1=(yy(ii+1,k1,i1,j1)-yy(ii,k1,i1,j1))/10.d0
width2=(pdf(i1,j1,ii+1)-pdf(i1,j1,ii))/10.d0
do jj=1,11
x(jj)=f0(ii)+width*(jj-1)
y1(jj)=yy(ii,k1,i1,j1)+width1*(jj-1)
p(jj)=pdf(i1,j1,ii)+width2*(jj-1)
fun(jj)=y1(jj)*p(jj)
enddo
call GenerateData(datas,width,func)
JF(jj)=Trape_Integral(fun,width)
sum=sum+JF(jj)
enddo
yym(k1,i1,j1)=sum
return
end subroutine