jixijifen.for

精细积分Fortran90程序

	IMPLICIT REAL*8 (A-H,O-Z)
	 real*8 h(4,4),f(4),hD(4,4),sum(4,4),ta(4,4),F0(4)
	data h/0.0d0,0.0d0,-6.0d0,2.0d0,0.0d0,0.0d0,2.0d0,-4.0d0,
	&       0.5d0,0.0d0,0.0d0,0.0d0,0.0d0,1.0d0,0.0d0,0.0d0/
	 data f/0.0d0,0.0d0,0.0d0,10.0d0/
	DATA F0/0.0d0,0.0d0,0.0d0,0.0d0/
	TAL=0.28D0
	DET=TAL/2**30
	HD=H*DET
	CALL REQURET(HD,DET,SUM,TA,4,30)
	ISTEP=0
19	 ISTEP=ISTEP+1
	F0=matmul(ta,f0)+matmul(sum,f) 
	WRITE(*,*)(F0(I),I=1,2)
	IF(ISTEP.LT.12)GOTO 19
	END
	SUBROUTINE REQURET(HD,DET,SUM,TA,N,NSTEP)
	IMPLICIT REAL*8 (A-H,O-Z)
	real(8),dimension(n,n):: HD,SUM,SUM1,TA
	ta=hd
	sum=hd/2
	sum1=hd
      do i=2,4
	sum1=matmul(sum1,hd) 
      ta=ta+sum1/i
	sum=sum+sum1/(i+1)
	enddo
	sum1=matmul(sum1,hd) 
	ta=ta+sum1/i
	 temp=2**nstep;i2=1
	DO 60 I=1,NSTEP
	sum=i2*ta+matmul(ta,sum)+2*sum
	ta=2*ta+matmul(ta,ta)
	i2=i2*2
60    continue
	DO 70 II=1,N
70	SUM(II,II)=SUM(II,II)+temp
	sum=sum*det
	DO 90 II=1,N
90    TA(II,II)=TA(II,II)+1.0d0
	END