## COPYRIGHT NOTICE ## Copyright (C) 2015 Edwin L. (Ted) Woollett ## http://www.csulb.edu/~woollett ## This program is free software; you can redistribute ## it and/or modify it under the terms of the ## GNU General Public License as published by ## the Free Software Foundation; either version 2 ## of the License, or (at your option) any later version. ## This program is distributed in the hope that it ## will be useful, but WITHOUT ANY WARRANTY; ## without even the implied warranty of MERCHANTABILITY ## or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details at ## http://www.gnu.org/copyleft/gpl.html ## FW.R uses Runge-Kutta for finite well. ## Ch. 3, Computational Physics with Maxima or R, ## Boundary Value and Eigenvalue Problems ## dimensionless units ## V = 1 for x < 0 and x > 1 ## V = 0 for 0 < x < 1 ## y''(x) + gam2*(E - V(x))*y(x) = 0 ## gam2 = gam^2 = 2500 ## gam = 50 = sqrt(2*m*L^2*V0/ hbar^2) ## initial global parameters: N = 100 h = 0.01 gam = 50 gam2 = gam^2 xdecay = 0.5 ## start yL1 integration at x = -xdecay ## start yR integration at x = 1 + xdecay ypleft = 1e-8 ypright = 1e-8 debug = FALSE wfdebug = FALSE cat (" gam = ",gam, " gam2 = ", gam2,"\n") cat (" h = ", h, ", xdecay = ", xdecay, ", ypleft = ", ypleft," ypright = ", ypright,"\n") ## analytic method functions Fa = function(k) { k - (2*k^2 - gam2)*tan(k)/2/sqrt(gam2 - k^2) } Frhs = function (k) { (2*k^2 - gam2)*tan(k)/2/sqrt(gam2 - k^2) } ## > Frhs(3) ## [1] 3.544391 ## start with fixed limits on plot to get the right idea: ## > curve (Frhs,1,gam - 1, n=200, col = "red",lwd = 3,xlim = c(0,gam),ylim = c(0,60), ## + xlab = "k", ylab = "") ## > lines( c(0,gam), c(0,gam),col = "blue",lwd = 3 ) kroot_plot = function (kkmin,kkmax, ymn, ymx) { curve (Frhs,kkmin,kkmax, n=200, col = "red",lwd = 3,ylim = c(ymn,ymx), xlab = "k", ylab = "") lines( c(kkmin,kkmax), c(kkmin,kkmax),col = "blue",lwd = 3 ) } kroot = function (kmin,kmax) uniroot (Fa, kmin, kmax, tol = 1e-10)$root ktoE = function (k) k^2/gam2 ## bracket_basic looks for a sign change in func, starting with x, ## and increasing x by dx each step. If sign change is found, ## then we back up to the previous x and search with new ## dx value one half of the previous value. ## Used by levels_analytic(kmax). ## debug set to FALSE at start of this file bracket_basic = function (func,xx,dxx,xacc) { x = xx dx = dxx it = 0 itmax = 1000 anerror = FALSE repeat { it = it + 1 if (it > itmax) { cat (" can't find change in sign \n") anerror = TRUE break} x1 = x x2 = x + dx if ( debug ) cat (" it = ",it," x1 = ",x1," x2 = ",x2," dx = ", dx, "\n") if (func(x1)*func(x2) < 0 ) { if ( abs(dx) < xacc ) break x = x - dx dx = dx/2 } else x = x2 } if (anerror) c(0,0) else c(x1,x2) } ## example using func = sin: ## > out = bracket_basic(sin,3,0.01,0.001) ## > out ## [1] 3.141250 3.141875 ## > uniroot(sin,out,tol=1e-10)$root ## [1] 3.141593 ## single root of x^3 ## > out = bracket_basic(function(x){x^3},-1,0.1,0.01) ## > uniroot(function(x) {x^3},out,tol=1e-10)$root ## [1] -1.387779e-16 ## zero node ground state using func = Fa: ## > out = bracket_basic(Fa,1.6,0.1,0.05) ## > out ## [1] 3.000 3.025 ## > kv = uniroot(Fa,out,tol = 1e-10)$root ## > kv ## [1] 3.02069 ## example can't find error ## > bracket_basic(function(x) {x^2},1,1,0.1) ## can't find change in sign ## [1] 0 0 ## analytic_wf(E) creates global functions yn1(x), yn2(x), yn0(x), yna(x) , ## and computes and prints analytic values of delx, delp. ## yn1(x) is the wave function for x < 0 ## yn2(x) is the wave function for x > 1 ## yn0(x) is the wave function for 0 <= x <= 1 ## yna(x) is the wave function for all x, usable by plot2d, but not by integrate. analytic_wf = function (E) { cat (" E = ", E, "\n") k = gam*sqrt(E) k1 = gam*sqrt( 1 - E ) k2 = k1 delta = atan(k/k1) sd = sin(delta) skd = sin(k + delta) D = sd^2/2/k1 + skd^2/2/k2 - sin(2*k+2*delta)/4/k + sin(2*delta)/4/k + 1/2 fac1 = 1/sqrt(D) A1 = fac1*sd A2 = fac1*skd ## calc < x > i1 = - A1^2/4/k1^2 iR = A2^2 * (1 + 2*k2)/4/k2^2 i2 = fac1^2*( -sin(2*k+2*delta)/(4*k)-cos(2*k+2*delta)/(8*k^2) + cos(2*delta)/(8*k^2)+1/4 ) x_mean = i1 + i2 + iR cat (" x_mean = ", x_mean, "\n" ) ## calc < x^2 > i1 = A1^2/4/k1^3 iR = A2^2 * (1 + 2*k2 + 2*k2^2)/4/k2^3 i2 = fac1^2 * ( -sin(2*k+2*delta)/(4*k)+sin(2*k+2*delta)/(8*k^3)-cos(2*k+2*delta)/(4*k^2) - sin(2*delta)/(8*k^3)+1/6 ) x2_mean = i1 + i2 + iR delx2 = x2_mean - x_mean^2 ## this should be positive! delx = sqrt(delx2) cat (" delx = ", delx, "\n" ) ## calc dyd_sum i1 = A1^2 / 2 iR = -A2^2 / 2 i2 = fac1^2 * ( cos(delta)^2 - cos(k + delta)^2 ) / 2 ydy_sum = i1 + i2 + iR cat (" ydy_sum = ", ydy_sum, "\n" ) ## calc i1 = -A1^2 * k1/2 iR = - A2^2 * k2/2 i2 = (fac1*k)^2 * ( - sin(2*k+2*delta)/4/k + sin(2*delta)/4/k + 1/2 ) delp2 = i1 + i2 + iR delp = sqrt(delp2) cat (" delp = ",delp," delx*delp = ", delx*delp, "\n" ) ## create globally usable functions ## yn1(x), yn2(x), yn0(x), yna(x) ## yna(x) usable with sapply(xL, yna) but not with curve(yna,0,1) yn1 <<- function (x) { A1*exp(k1*x) } yn2 <<- function (x) {A2*exp(k2)*exp(-k2*x) } yn0 <<- function (x) { fac1*sin(k*x + delta) } yna <<- function (x) { if (x < 0) yn1(x) else if (x > 1) yn2(x) else yn0(x) } } ## ## > kv ## [1] 3.02069 ## > Ev = ktoE(kv) ## > Ev ## [1] 0.00364983 ## > analytic_wf(Ev) ## E = 0.00364983 ## x_mean = 0.5 ## delx = 0.18802 ## ydy_sum = -1.78677e-16 ## delp = 2.96193 delx*delp = 0.556902 ## > xL = seq(-0.5,1.6,by = 0.01) ## > yL = sapply(xL, yna) ## > fll(yL) ## 1.21784e-12 8.28101e-15 211 ## > fll(xL) ## -0.5 1.6 211 ## > nodes_analytic(0.01) ## [1] 0 ## > plot(xL,yL,type = "l",lwd = 3, col = "blue",xlab = "x",ylab = "y(x)",tck=1) ## > plot_analytic(Ev) ## E = 0.00364983 ## x_mean = 0.5 ## delx = 0.18802 ## ydy_sum = -1.78677e-16 ## delp = 2.96193 delx*delp = 0.556902 ## ymn = 0 ymx = 2 ## number of nodes = 0 ## nodes_analytic(dx) counts the number of nodes implied by ## the function yn0(x) created by analytic_wf(E) nodes_analytic = function ( dx ) { num = 0 xv = 0 repeat { f1 = yn0(xv) xnew = xv + dx if (xnew > 1) break f2 = yn0(xv + dx) if (f1*f2 < 0 ) num = num + 1 xv = xnew } num } ## plot_analytic(E) calls analytic_wf(E), prints out ## the number of nodes, and makes a plot of yna(x) plot_analytic = function (E) { ddx = 0.001 xmn = - 0.25 xmx = 1.25 analytic_wf(E) xL = seq(xmn,xmx ,by = ddx) yL = sapply(xL, yna) ymn = floor( min(yL) ) ymx = 1 + floor( max (yL) ) cat (" ymn = ",ymn," ymx = ",ymx, "\n") cat (" number of nodes = ", nodes_analytic(ddx), "\n" ) plot(xL, yL, type = "l", lwd = 3, col = "blue", xlab = "x",ylab = "ya(x)",tck=1, ylim = c(ymn,ymx) ) } ## > plot_analytic(Ev) ## E = 0.00364983 ## x_mean = 0.5 ## delx = 0.18802 ## ydy_sum = -1.78677e-16 ## delp = 2.96193 delx*delp = 0.556902 ## ymn = 0 ymx = 2 ## number of nodes = 0 ## levels_analytic(kmax) returns a vector of analytic eigen energies ## related to k via: k = gam*sqrt(E). The corresponding eigenvalues ## of k are separated by roughly dk = 3 and lie in the middle of the intervals ## [n1*%pi/2, n2*%pi/2], where n1 and n2 are adjacent odd integers, n2 > n1. ## The ground state (zero node soln) corresponds to k = 3.02 (n1=1,n2=3). levels_analytic = function (kmax) { rmax = 20 EL = rep(NA, rmax) level = 0 nmx = ceiling(2*kmax/pi) ## make nmx an odd integer if ( is.even (nmx) ) nmx = nmx + 1 ## cat (" nmx = ", nmx, "\n") cat (" nodes E \n ") cat ( " \n ") ## analytic kv using n*pi/2 + 0.1 with n odd for ( j in seq(1, nmx, by=2) ) { out = bracket_basic( Fa, j*pi/2 + 0.1, 0.01,0.005) kv = uniroot( Fa, out, tol = 1e-16)$root Ev = ktoE(kv) EL[ j ] = Ev cat ( " ", level, " ", Ev, "\n" ) level = level + 1 } EL[!is.na(EL)] } ## > levels_analytic(20) ## nodes E ## 0 0.00364983 ## 1 0.0145973 ## 2 0.032836 ## 3 0.0583552 ## 4 0.0911389 ## 5 0.131165 ## 6 0.178405 ## [1] 0.00364983 0.01459726 0.03283598 0.05835517 0.09113890 0.13116525 0.17840502 ### functions for numerical integration ## wf(E) creates ** un-normalized ** numerical wave functions ## using Runge-Kutta routine myrk4. ## The wave functions are stored in global vectors ## xL1, yL1,ypL1, xL2, yL2, ypL2, xR, yR, ypR . ## Program also defines **global** nleft, nright, ncenter ## the global xL1 grid extends from -xdecay to 0 and ## the global xL2 grid extends from 0 to 1 and ## the global xR grid extends from 1 to 1 + xdecay wf = function (E) { if ( (E < 0) | (E > 1) ) { cat (" need 0 < E < 1 \n") return(false) } ncenter = N if ( round(ncenter) != ncenter) { cat (" ncenter = ",ncenter," is not an integer \n") return(false) } nleft = round(xdecay/h) nright = nleft if (wfdebug) cat (" nleft = ",nleft," ncenter = ",ncenter," nright = ",nright, "\n") glr = gam2*(E - 1) ## g(x) for x < 0 and x > 1 gc = gam2*E ## g(x) for 0 < x < 1 if (wfdebug) cat (" glr = ",glr," gc = ", gc, "\n") ## construct xl1, yl1, and ypl1 for -xdecay < x < 0 derivs.decay = function (x,y) { c (y[2], - glr*y[1]) } xl1 = seq (- xdecay, 0, h) outL = myrk4( c(0,ypleft), xl1, derivs.decay) yl1 = outL[[1]] ypl1 = outL[[2]] ## construct xl2, yl2, and ypl2 for 0 < x < 1 derivs01 = function (x,y) { c (y[2], - gc*y[1]) } xl2 = seq(0, 1, h) outL = myrk4( c(last(yl1), last(ypl1)), xl2, derivs01) yl2 = outL[[1]] ypl2 = outL[[2]] ## construct xr, yr, and ypr for 1 < x < 1 + xdecay xr = seq( 1 + xdecay, 1, -h) outL = myrk4( c(0, -ypright), xr, derivs.decay) yr = outL[[1]] ypr = outL[[2]] xr = rev(xr) yr = rev(yr) ypr = rev(ypr) if (wfdebug) cat (" yr(1) = ", yr[1], "\n" ) fac = yr[1] / last(yl2) if (wfdebug) cat (" fac = ",fac, "\n") yl1 = fac*yl1 ypl1 = fac*ypl1 yl2 = fac*yl2 ypl2 = fac*ypl2 ## make global xL1,xL2,xR,yL1,yL2,yR,ypL1,ypL2,ypR xL1 <<- xl1 xL2 <<- xl2 xR <<- xr yL1 <<- yl1 yL2 <<- yl2 yR <<- yr ypL1 <<- ypl1 ypL2 <<- ypl2 ypR <<- ypr } ## > source("FW.R") ## gam = 50 gam2 = 2500 ## h = 0.01 , xdecay = 0.5 , ypleft = 1e-08 ypright = 1e-08 ## > wf(0.5) ## > fll(xL1) ## -0.5 0 51 ## > fll(yL1) ## 0 -0.00475513 51 ## > fll(ypL1) ## -7.08074e-09 -0.168119 51 ## > fll(xL2) ## 0 1 101 ## > fll(yL2) ## -0.00475513 0.00671558 101 ## > fll(ypL2) ## -0.168119 -0.00190471 101 ## > fll(xR) ## 1 1.5 51 ## > fll(yR) ## 0.00671558 0 51 ## > fll(ypR) ## -0.237432 -1e-08 51 ## dy_diff() uses global ypL2, ypR, and yR, ## returns a normalized difference of derivatives ## (yL'(1) - yR'(1)/ yR(1) dy_diff = function () { dy_left = last(ypL2) dy_right = ypR[1] (dy_left - dy_right)/abs( yR[1] ) } ## energy eigenvalue if global function F(E) = 0 . ## F(E) calls wf(E) then returns dy_diff(), but ## returns FALSE if E < 0 or > 1 . F = function (E) { if (E < 0 | E > 1) { cat (" in F(E), E = ",E," should be between 0 and 1 \n ") return(FALSE) } wf(E) dy_diff() } ## > F(0.5) ## [1] 35.0717 ## > F(1.2) ## in F(E), E = 1.2 should be between 0 and 1 ## [1] FALSE ## bracket is a modified version of bracket_basic, designed to work with ## the functions F(E) or F1(k) which can return FALSE. ## bracket looks for a sign change in func, ## starting with xx, and increasing xx by dxx each step. ## If sign change is found, then we back up to the previous xx ## and search with new dxx value one half of the previous value. ## normally returns [ea,eb], but if can't find change in sign, ## then returns [0,0], and if func returns FALSE, then ## bracket returns FALSE. bracket = function (func,xx,dxx,xacc) { x = xx dx = dxx it = 0 itmax = 1000 anerror = FALSE anerror2 = FALSE repeat { it = it + 1 if (it > itmax) { cat (" can't find change in sign \n") anerror = TRUE break} x1 = x x2 = x + dx if ( debug ) cat (" it = ",it," x1 = ",x1," x2 = ",x2," dx = ", dx, "\n") f1 = func(x1) if ( f1 == FALSE) { cat (" in bracket, f1 = FALSE , x1 = ",x1, " dx = ", dx, " \n ") anerror2 = TRUE break } f2 = func(x2) if ( f2 == FALSE) { cat (" in bracket, f2 = FALSE , x2 = ",x2, " dx = ", dx, " \n ") anerror2 = TRUE break } if ( f1 * f2 < 0 ) { if ( abs(dx) < xacc ) break x = x - dx dx = dx/2 } else x = x2 } if (anerror) c(0,0) else if (anerror2) FALSE else c(x1,x2) } ## > out = bracket(F,0.0005,0.0001,0.00005) ## > out ## [1] 0.003625 0.003650 ## > e = uniroot(F,out, tol = 1e-10)$root ## > e ## [1] 0.00364983 ## count the number of nodes in yL2 num_nodes = function () { n = 0 for ( j in 1 : (length(yL2) - 1) ) { if (yL2[ j ] * yL2[ j + 1 ] < 0) n = n + 1 } n } ## normalize() uses the current global xL1,yL1, xL2,yL2, xR, yR and ## the utility functions rest and simp to define global ## xn and yn, with the latter being normalized. normalize = function () { xn = c (xL1, c ( rest(xL2), rest(xR) ) ) yn = c (yL1, c ( rest(yL2), rest(yR) ) ) ## we need xn to have odd number of elements to use simp if ( is.even ( length (xn) ) ) { xn = rest (xn) yn = rest (yn) } if (debug) cat ( " fll(xn) = ", fll(xn), "\n" ) if (debug) cat ( " fll(yn) = ", fll(yn), "\n" ) AA = simp(xn,yn^2) cat ( " AA = ",AA, "\n" ) yn = yn/sqrt(AA) if (debug) cat ( " fll(yn) = ", fll(yn), "\n" ) x_mean = simp(xn, xn * yn^2) cat (" x_mean = ", x_mean, "\n" ) x2_mean = simp(xn, xn^2 * yn^2) delx2 = x2_mean - x_mean^2 ## this should be positive! delx = sqrt(delx2) cat (" delx = ", delx, "\n" ) xn <<- xn yn <<- yn } ## yn_plot_current() uses the currently defined normalized set (xn,yn) yn_plot_current = function () { ymn = floor( min(yn) ) ymx = 1 + floor( max (yn) ) cat (" ymax = ", max(yn), "\n" ) plot(xn, yn, type = "l", lwd = 3, col = "blue", ylim = c(ymn, ymx), xlab = "x", ylab = "y", tck = 1) } ## yn_plot(E,xmin,xmax) first calls wf(E) to create ## un-normalized wave functions corresponding to the ## given energy E. Then normalizes those wave functions ## to produce the vectors xn and yn. Finally makes a plot ## of yn over only the region (xmn, xmx) yn_plot = function (E,xmn,xmx) { wf(E) cat (" E = ",E, "\n" ) cat (" number of nodes = ",num_nodes(),", dy_diff = ",dy_diff(), "\n" ) normalize() cat (" normalized ymax = ", max(yn), "\n" ) ymn = floor( min(yn) ) ymx = 1 + floor( max (yn) ) plot(xn, yn, type = "l", lwd = 3, col = "blue", xlim = c(xmn,xmx), ylim = c(ymn, ymx), xlab = "x", ylab = "y", tck = 1) } ## F1(k) ## energy eigenvalue if global function F1(k) = 0 . ## F1(k) calls wf(E) then returns dy_diff(), but ## returns false if k <= 0 or k >= gam. F1 = function (k) { if (k <= 0 | k >= gam) { cat (" in F1(k), k = ",k," k should be greater than 0 and less than ",gam, "\n") return(FALSE) } wf(k^2/gam2) dy_diff() } ## > F1(gam+1) ## in F1(k), k = 51 k should be greater than 0 and less than 50 ## [1] FALSE ## > F1(-1) ## in F1(k), k = -1 k should be greater than 0 and less than 50 ## [1] FALSE ## levels(kmin,kmax,dk, kacc ) returns a vector c( Ea, Eb,...) of energy levels with ## increasingly larger number of nodes in energy range (Emin, Emax) ## according to Emin = kmin^2/gam^2, and Emax = kmax^2/gam^2. ## uses F1(k) (inside bracket) to find roots, and calls wf(E), num_nodes() and dy_diff() for each root found, ## Uses bracket and uniroot. ## The arguments (dk, kacc) are used to call bracket, and do not describe ## the accuracy of the energy levels found. ## Once a good energy e.v. is found we look for the region of ## energies with one more node and search there. ## Code has interactive continue or stop. ## Searching for the k eigenvalues via F1(k) is easier than searching ## directly for the E eigenvalues via F(E) for the case gam = 50 we ## consider in this example. levels = function (kmin,kmax,dk, kacc ) { rmax = 20 eL = rep(NA, rmax) ## vector eL will hold energy eigenvalues found k = kmin nlast = -1 j = 1 repeat { if (k > kmax | j > rmax) break ## exit do loop cat ("---------------- levels -------------------\n") cat (" nlast = ", nlast,"\n") cat (" kstart = ", k," dk = ", dk, "\n" ) out = bracket(F1,k, dk, kacc) cat (" ka = ",out[1]," kb = ",out[2],"\n") if (out[1] == 0) { cat (" can't find bracket interval \n") cat (" k = ",k, "\n") break } kroot = uniroot(F1, out, tol = 1e-16)$root cat (" kroot = ", kroot, "\n") eroot = kroot^2/gam2 cat (" eroot = ", eroot, "\n") wf(eroot) nn = num_nodes() cat (" number of nodes = ", nn, "\n") cat (" dy_diff at x = 1 is ", dy_diff(), "\n" ) eL[ j ] = eroot nlast = nn j = j + 1 r = readline (" input c or s \n ") if (r == "s") break ## exit do loop ## search for a k value greater than kb which produces ## a wave function with nn + 1 nodes knext = out[2] + dk repeat { wf(knext^2/gam2) if (num_nodes() > nlast) { k = knext break } else knext = knext + dk } } ## end of outer repeat loop ## remove na's at end of eL eL[!is.na(eL)] } tryit = function() { j = 1 repeat { if (j > 5) break r = readline(" input c or s \n ") if (r == "s") { cat( " r = s \n") break } cat ( " r = c \n") j = j + 1 } cat(" bye \n") } ## > tryit() ## input c or s ## c ## r = c ## input c or s ## c ## r = c ## input c or s ## s ## r = s ## bye