Re: Numerical Solutions to Schrod Eq with Mathematica
- To: mathgroup at smc.vnet.net
- Subject: [mg2365] Re: Numerical Solutions to Schrod Eq with Mathematica
- From: lamontg at u.washington.edu (Lamont Granquist)
- Date: Fri, 27 Oct 1995 01:51:16 -0400
- Organization: University of Washington
lamontg at u.washington.edu (Lamont Granquist) writes:
>Okay, so I've got Mathematica plugging away at the Numerov method of
>solving the time-independent schrodinger equation for even potentials.
>I can get the lower eigenvalues. Now, the trick is to plug them into
>NDSolve with the appropriate initial conditions. I've been only playing
>with even solutions to SE for now, so I can pick psi'[0] == 0 as one
>BC, but then the best I can do is to approximate psi[0] by the second
>term in the appropriate eigenvector (for some reason the first term in
>the eigenvector always seems to be corrupt). Is there a better way to
>get psi[0]? Is there a better choice of initial conditions?
Okay, so *duh*, the psi[0] value is nothing other than the normalization
constant. It probably works better to get the normalization constant
right on before handing it to NDSolve, but it shouldn't cause that much
error unless it is grossly off. It also looks sucpiciously a lot like
the first term in the eigenvector is just 1/2 of psi[0], which makes
a lot of sense.
For anyone who might be interested in this kind of thing here's the
mathematica code I've got so far on this problem. One thing I'd really
like some help on is finding out how to generate the terms for the
expansion of D^2 so that I could write a program which would automatically
generate the nth extended numerov method.
Needs["LinearAlgebra`MatrixManipulation`"]
EvenAMatrix[P_, h_, R_] := Module[
{AMatrix, N, V} ,
N = Floor[R/h] ;
V[k_] = P /. x -> k h ;
AMatrix = ZeroMatrix[N, N] ;
AMatrix[[1,1]] = -24 - 10 h^2 V[0] ;
AMatrix[[2,1]] = 2 ( 12 - h^2 V[1] ) ;
AMatrix[[N-1,N]] = 12 - h^2 V[N-2] ;
AMatrix[[N,N]] = -24 - 10 h^2 V[N-1] ;
For[i = 1, i < N-1, i = i + 1,
AMatrix[[i,i+1]] = 12 - h^2 V[i-1] ;
AMatrix[[i+1,i+1]] = -24 - 10 h^2 V[i] ;
AMatrix[[i+2,i+1]] = 12 - h^2 V[i+1] ;
] ;
AMatrix
] ;
EvenBMatrix[P_, h_, R_] := Module[
{BMatrix, N, V} ,
N = Floor[R/h] ;
V[k_] = P /. x -> k h ;
BMatrix = ZeroMatrix[N, N] ;
BMatrix[[1,1]] = 10 ;
BMatrix[[2,1]] = 2 ;
BMatrix[[N-1,N]] = 1 ;
BMatrix[[N,N]] = 10 ;
For[i = 1, i < N-1, i = i + 1,
BMatrix[[i,i+1]] = 1 ;
BMatrix[[i+1,i+1]] = 10 ;
BMatrix[[i+2,i+1]] = 1 ;
] ;
BMatrix
] ;
A = EvenAMatrix[ x^2, .02, 4] //N ;
Print["found matrix A"]
B = EvenBMatrix[ x^2, .02, 4] //N ;
Print["found matrix B"]
BInv = Inverse[B] //N ;
Print["found inverse of B"]
BInvA = BInv . A //N ;
Print["found B^-1A"]
Evect = Eigensystem[BInvA] //N ;
Print["found eigensystem"]
--
Lamont Granquist (lamontg at u.washington.edu)
There comes a point, I'm afraid, where you begin to suspect that if there's
any _real_ truth, it's that the entire multidimentional infinity of the
Universe is almost certainly being run by a bunch of maniacs -- Douglas Adams