메뉴 건너뛰기

나는 갈 데가 없다

삶

IT 컴공이나 수학과 있으면 Fortran 에러좀 잡아줘

2014.06.19 17:20

리카아메 조회 수:670

program oscillator
implicit none
integer :: L,M
integer :: i,j,ell
real(8),dimension(:),allocatable :: u,vpot,Ei
real(8) :: h,x_b,kappa,E,gamma_4
real(8) :: x,tmp_E
real(8) :: c1,c2
real(8) :: f0,f1,f2
real(8) :: El,Eu
L=8192
M=100
kappa=1.0D0
gamma_4=1.0D0
x_b=dsqrt(36.8/dsqrt(kappa))
h=x_b/dble(L)
allocate(u(0:L-1))
allocate(vpot(0:L-1))
allocate(Ei(0:M))
! ---initial guess for lower and upper bound of E---
El=0.0D0
Eu=5.0D0
      Ei(0)=0.5d0*(El+Eu)

c1=h*h/12.0D0
c2=5.0D0*h*h/6.0D0

! ---iteration---
do i=0,M-1

tmp_E=Ei(i)
! ---initial guess for a ground sate---
u(0)=1.0D0
u(1)=1.0D0-0.5D0*tmp_E*(h**2)
!u(0)=0.0D0
!u(1)=1.0D0*h

! ---potential---
do ell=0,L-1

x=dble(ell)*h
vpot(ell)=(kappa+gamma_4*x*x)*x*x-tmp_E

enddo

! ---Numerov---
do ell=2,L-1

f2 = 1.0D0 - c1*vpot(ell-2)
f1 = 2.0D0 + c2*vpot(ell-1)
f0 = 1.0D0 - c1*vpot(ell)
u(ell)=(f1*u(ell-1) - f2*u(ell-2))/f0

enddo

write(11,*)i,Ei(i)

! ---update El and Eu
if(u(L-1).gt.1.0D-5) then
El=Ei(i)
Ei(i+1)=0.5D0*(El+Eu)
else if(u(L-1).lt.(-1.0D-5)) then
Eu=Ei(i)
Ei(i+1)=0.5D0*(El+Eu)
else
write(*,*)"# converge!",i
exit
endif

enddo

deallocate(u,vpot,Ei)
end program oscillator



시발 아무리 해도 if에서 제대로 converge했다고 안뜨는데 도대체 왜이런거임?
결과 봐도 u(L-1)은 진즉에 수렴했는데 도대체 왜...


과제라서 급한데 연구실에 다 학회가서 물어볼인간이 없네;;




이 게시물을
<?php if(!defined("__XE__")) exit();?>

                                                                 영원불멸의 프리큐어5


R2VY3L2.png?1

목록
목록
쓰기
태그
위로