c Problem 3 for colnew, example of continuation
c Quantum Neumann equation, with 2 "eigenvalues" c_1 and c_2
c (c_1=v-c_2-c_3, v is a parameter, used in continuation)
c
c diff(f,x,2) + (1/2)*(1/x + 1/(x-1) + 1/(x-y))*diff(f,x)
c - (c_1/x + c_2/(x-1) + c_3/(x-y))* f(x) = 0
c
c with 4 "boundary" conditions: diff(f,x)(a_k)=2*c_k*f(a_k) for
c k=1,2,3, a_k=(0, 1, y) and normalization f(1) = 1
c
c The z-vector is z_1=f, z_2=diff(f,x), z_3=c_2 and z_4=c_3
c The guess is chosen to have one node in [0,1], f(x)=2*x-1
c such that f(1)=1, c_2 and c_3 are chosen to cancel poles in
c the differential equation at 1.0 and y, z_3=1, z_4=1/(2*y-1)
subroutine fsub(x, z, f)
common /cprob3/v, y
real*8 v, y
real*8 z(4), f(3), x
f(1) = - (1/2.0d0)*(1/x + 1/(x-1.0d0) + 1/(x-y))*z(2) +
$ ((v - z(3) - z(4))/x + z(3)/(x-1.0d0) + z(4)/(x-y))*z(1)
f(2) = 0
f(3) = 0
return
end
subroutine dfsub(x, z, df)
common /cprob3/v, y
real*8 v, y
real*8 z(4), f(3), x, df(3, 4)
df(1,1) = z(4)/(x-y)+(-z(4)-z(3)+v)/x+z(3)/(x-1.0d0)
df(1,2) =-0.5d0*(1/x+1/(x-1.0d0)
$ +1/(x-1.9d0))
df(1,3) = (1/(x-1.0d0)-1/x)*z(1)
df(1,4) =(1/(x-y)-1/x)*z(1)
df(2,1) =0
df(2,2) =0
df(2,3) =0
df(2,4) =0
df(3,1) =0
df(3,2) =0
df(3,3) =0
df(3,4) =0
c print *, ' df =', df
return
end
subroutine gsub(i, z, g)
common /cprob3/v, y
real*8 v, y
real*8 z(4), g
go to (10, 20, 30, 40), i
10 g = z(2) - 2.0d0*z(1)*(v-z(3)-z(4))
go to 99
20 g = z(2) - 2.0d0*z(1)*z(3)
go to 99
30 g = z(1) - 1.0d0
go to 99
40 g = z(2) - 2.0d0*z(1)*z(4)
99 continue
return
end
subroutine dgsub(i, z, dg)
common /cprob3/v, y
real*8 v, y
real*8 z(4), dg(4)
go to (10, 20, 30, 40), i
10 dg(1) = -2.0d0*(-z(4)-z(3)+v)
dg(2) = 1.0d0
dg(3) = 2.0d0*z(1)
dg(4) = 2.0d0*z(1)
go to 99
20 dg(1) = -2.0d0*z(3)
dg(2) = 1.0d0
dg(3) = -2.0d0*z(1)
dg(4) = 0
go to 99
30 dg(1) = 1.0d0
dg(2) = 0
dg(3) = 0
dg(4) = 0
go to 99
40 dg(1) = -2.0d0*z(4)
dg(2) = 1.0d0
dg(3) = 0
dg(4) = -2.0d0*z(1)
99 continue
return
end
subroutine guess(x, z, dmval)
common /cprob3/v, y
real*8 v, y
real*8 x, z(4), dmval(3)
z(1) = 2.0d0*x-1.0d0
z(2) = 2.0d0
z(3) = 1.0d0
z(4) = 1/(2.0d0*y-1.0d0)
dmval(1) = 0
dmval(2) = 0
dmval(3) = 0
return
end
program prob3
implicit real*8 (a-h, o-z)
real*8 fspace(30000), zeta(4), tol(4), z(4)
integer ispace(2000), m(3), ipar(11), ltol(4)
external fsub, dfsub, gsub, dgsub, guess
real*8 fixpnt(5)
common /cprob3/v, y
real*8 v, y
parameter (nsteps=100)
real*8 eig(nsteps, 4)
y = 1.9d0
v = 0.00001d0
m(1) = 2
m(2) = 1
m(3) = 1
zeta(1) = 0.0d0
zeta(2) = 1.0d0
zeta(3) = 1.0d0
zeta(4) = y
do 10 i = 1, 11
ipar(i) = 0
10 continue
c Non-linear prob
ipar(1) = 1
c Number of collocation points
ipar(2) = 3
c Initial uniform mesh of 4 subintervals
ipar(3) = 4
ipar(8) = 0
c Size of fspace, ispace
ipar(5) = 30000
ipar(6) = 2000
c medium output, set it to 1 to have a clean table of eigenvalues
ipar(7) = 0
c Initial approx provided
ipar(9) = 1
c fixpnt is an array containing all fixed points in the mesh
ipar(11) = 1
c tolerances on all components
ipar(4) = 4
c tolerance on f and diff(f,x) and on c2 and c3
do 20 i = 1, 4
ltol(i) = i
tol(i) = 1d-5
20 continue
c only one interior boundary
fixpnt(1) = 1.0d0
c First run, using the guess
call colnew(3, m, 0d0, y, zeta, ipar, ltol, tol, fixpnt, ispace,
$ fspace, iflag, fsub, dfsub, gsub, dgsub, guess)
len = int(1+y/.1d0)
print *, 'Wave function at v=0'
do 30 i = 0, len
call appsln(.1d0*i, z, fspace, ispace)
c Print z(1) to see the solution here, a straight line
print *, z(1)
30 continue
print *, 'Eigenvalues'
write(*, 9999) v, z(3), z(4), z(3)+y*z(4)
c Perform continuation, which implies modifying ipar(9) and ipar(3)
ipar(9) = 3
ipar(8) = 1
do 100 i = 1, nsteps
ipar(3) = ispace(1)
c poor man logarithmic scale
if ( i .lt. nsteps/3) v = v + 0.2d0
if ( i .ge. nsteps/3) v = v + 0.4d0
if ( i .gt. 2*nsteps/3) v = v + 0.6d0
call colnew(3, m, 0d0, y, zeta, ipar, ltol, tol, fixpnt,
$ ispace,
$ fspace, iflag, fsub, dfsub, gsub, dgsub, guess)
call appsln(.1d0, z, fspace, ispace)
eig(i,1) = z(3)
eig(i,2) = z(4)
eig(i,3) = z(3)+y*z(4)
eig(i,4) = v
100 continue
do 110 i = 1, nsteps
write(*,9999) eig(i, 4), eig(i, 1), eig(i, 2), eig(i, 3)
110 continue
9999 format(5H v = ,f6.3,7H c_2 = ,f6.3,7H c_3 = ,f6.3,7H E/2 = ,f6.3)
c Show the solution for the last value of v
print *, 'Wave function at large v'
do 120 i = 0, len
call appsln(.1d0*i, z, fspace, ispace)
print *, z(1)
120 continue
c The results show that the probability of presence is concentrated
c near t=1 and near t=y which corresponds to a pole of the sphere.
end