Solving a QM problem numerically

[rahuldandekar]

I've been trying to solve the Deuteron Tensor polarisation equations on the PC, in Python, using a fourth-order Runge-Kutta method. I have a couple of questions:

1) Is Python ok for this job, in that does it give the right answers? The speed is not an issue, it's fast enough for me (I use Numpy and psyco).

2) (which is much more important) How do I cancel out the solutions that "blow up" at infinity? When the computer solves the equations, it gives the most general solution, and doesn't care that the solution blows up at infinity. I DO care ;) , and would like to only get the solution which doesn't blow up... and that would give me th right answer... how do I do it?

 

[jim mcnamara]

1. Python uses double precision floating point to represent numbers. Only you can decide if that meets your needs. Double precision gives 15 signficant digits in a result.
See:
http://steve.hollasch.net/cgindex/coding/ieeefloat.html

2. I don't understand what you mean by cancel out -discard the result?

 

[George Jones]

2) (which is much more important) How do I cancel out the solutions that "blow up" at infinity? When the computer solves the equations, it gives the most general solution, and doesn't care that the solution blows up at infinity. I DO care ;) , and would like to only get the solution which doesn't blow up... and that would give me th right answer... how do I do it?

A numerical solution to a (system of) differential equation(s), is not a general solution, it is (a close approximation to) the *unique* solution for the initial values that you specify.

If your numerical solution is blowing up, possibilities include:

1) you have specified the wrong initial conditions;

2) you haven't coded things correctly;

3) the solution depends critically on initial values, with very small variations in initial conditions producing wildly different solutions.

 

[rahuldandekar]

Thanks for the help. :) In fact double precision is what the book says is best... it allows one to use smaller intervals.

About blowing up, what I mean is this:

Usually, you get two solutions to an equation... say e^x + e^-x. Now, in QM, I want the solution to be normalisable, and thus I discard the e^x solution because it doesn't go to 0 at +infinity. However, when you solve numerically, the computer will see the initial conditions, and give an appropriate complete solution, NOT discarding the e^x term. Thus the numerical solution is of no use to me if it's not normalisable.

There must be a way around this... I've tried some trickery, but it doesn't work.

 

[George Jones]

I don't quite see what you're after, so let me try out a toy example.

The differential equation

y'' - y = 0

has general solution

y(x) = a e^x + b e^{-x}.

Suppose the initial conditions are y(0) = 1 and y'(0) = 0.

This implies a = b = 1/2, so, for these initial conditions,

y(x) = 1/2 (e^x + e^{-x}).

Throwing away the e^x term as unphysical gives

y(x) = 1/2 e^{-x},

but this isn't a solution of the given differential equation for the given initial conditions.

For this to be a solution, either the differential equation or the initial conditions, or both, has to change.

For example, keeping the original DE, but changing the initial conditions to y(0) = 1/2 and y'(0) = -1/2 results in the solution

y(x) = 1/2 e^{-x}.

Maybe I'm way off base, but I think choosing appropriate initial conditions might eliminate exponentially growing solutions.

 

[mda]

Maybe I'm way off base, but I think choosing appropriate initial conditions might eliminate exponentially growing solutions.

I tend to agree that the boundary conditions are important here.

I'm not familiar with the particular system under study, but I'd suggest that systems in the spatial domain (including QM) should be initially cast as BVP not IVP.
In practice this means approximating the normalisation condition by explicitly requiring the function at large distance to be zero. As a starting point investigate the 'shooting method'.