CODEWITHSUNDEEP
pythonsympysymbolic-mathequation-solving
Galen
Galen

Reputation: 1352

Why won't SymPy directly solve for the normalization of this wavefunction?

Here is some minimal code demonstrating the problem.

from sympy import *

x = Symbol('x', real=True)
t = Symbol('t', real=True)
A = Symbol('A', real=True, positive=True)
λ = Symbol('λ', real=True, positive=True)
ω = Symbol('ω', real=True, positive=True)

# Define wavefunction
psi_x_t = A * exp(-λ * Abs(x)) * exp(-I*ω*t)

# Normalize
print(solve(integrate(psi_x_t**2, (x, -oo, oo)) - 1))
print(solve(integrate(psi_x_t**2, (x, -oo, oo)) - 1, A))

Where I simply call solve it returns [{t: 0, λ: A**2}]. What I am looking for is actually A=sqrt(λ). But when I include A as the symbol I am solving for I don't get a plus/minus sqrt(λ) but rather an empty collection [].

Why is solve giving me an empty list of solutions, and is there a way to make it directly return the solution for A?

Upvotes: 2

Views: 97

Answers (3)

Galen
Galen

Reputation: 1352

It did not solve for A as expected because it doesn't choose some time in which to solve for it. If we first set time, such as t=0, then it will find a suitable solution.

from sympy import *

x = Symbol('x', real=True)
t = Symbol('t', real=True)
A = Symbol('A', real=True, positive=True)
λ = Symbol('λ', real=True, positive=True)
ω = Symbol('ω', real=True, positive=True)

# Define wavefunction
psi_x_t = A * exp(-λ *Abs(x)) * exp(-I*ω*t)

# Integrate
result = integrate(psi_x_t**2, (x, -oo, oo)) - 1
result = result.subs(t, 0)
# Normalize
print(solve(result))
print(solve(result, A))

Upvotes: 0

Oscar Benjamin
Oscar Benjamin

Reputation: 14500

You are not getting a solution because you have declared assumptions on your variables and those assumptions are inconsistent with the solution. Let's try without those assumptions:

In [20]: from sympy import *
    ...: 
    ...: x = Symbol('x')
    ...: t = Symbol('t')
    ...: A = Symbol('A')
    ...: λ = Symbol('λ')
    ...: ω = Symbol('ω')
    ...: 
    ...: # Define wavefunction
    ...: psi_x_t = A * exp(-λ * Abs(x)) * exp(-I*ω*t)

In [21]: psi_x_t
Out[21]: 
   -λ⋅│x│  -ⅈ⋅t⋅ω
A⋅ℯ      ⋅ℯ      

In [22]: eq = integrate(psi_x_t**2, (x, -oo, oo)) - 1

In [23]: eq
Out[23]: 
⎛⎧         2  -2⋅ⅈ⋅t⋅ω                          ⎞    
⎜⎪        A ⋅ℯ                                 π⎟    
⎜⎪        ────────────          for │arg(λ)│ < ─⎟    
⎜⎪             λ                               2⎟    
⎜⎪                                              ⎟    
⎜⎨∞                                             ⎟ - 1
⎜⎪⌠                                             ⎟    
⎜⎪⎮   2  -2⋅λ⋅│x│  -2⋅ⅈ⋅t⋅ω                     ⎟    
⎜⎪⎮  A ⋅ℯ        ⋅ℯ         dx     otherwise    ⎟    
⎜⎪⌡                                             ⎟    
⎝⎩-∞                                            ⎠

I'm guessing that this Piecewise is the reason that you set those assumptions. Let's just manually select the first case of the Piecewise and work from there:

In [26]: eq1 = piecewise_fold(eq).args[0][0]

In [27]: eq1
Out[27]: 
 2  -2⋅ⅈ⋅t⋅ω    
A ⋅ℯ            
──────────── - 1
     λ

Now we solve this:

In [28]: solve([eq1], [A])
Out[28]: 
⎡⎛     ⅈ⋅t⋅ω ⎞  ⎛    ⅈ⋅t⋅ω ⎞⎤
⎣⎝-√λ⋅ℯ     ,⎠, ⎝√λ⋅ℯ     ,⎠⎦

Now for almost all possible values of omega this is not going to be real which is a violation of the assumptions you stated that A should be positive.

Presumably the equation you actually wanted comes from integrating abs(psi)^2:

In [34]: from sympy import *
    ...: 
    ...: x = Symbol('x', real=True)
    ...: t = Symbol('t', real=True)
    ...: A = Symbol('A', real=True, positive=True)
    ...: λ = Symbol('λ', real=True, positive=True)
    ...: ω = Symbol('ω', real=True, positive=True)
    ...: 
    ...: # Define wavefunction
    ...: psi_x_t = A * exp(-λ * Abs(x)) * exp(-I*ω*t)

In [35]: eq = integrate(abs(psi_x_t)**2, (x, -oo, oo)) - 1

In [36]: eq
Out[36]: 
 2    
A     
── - 1
λ     

In [37]: solve([eq], [A])
Out[37]: [(√λ,)]

Here only the positive root is returned because we assumed A to be positive.

Upvotes: 2

Davide_sd
Davide_sd

Reputation: 13185

Why is solve giving me an empty list of solutions, and is there a way to make it directly return the solution for A?

I can't answer this question. However, when solve fails I usually give a try to solveset:

print(solveset(integrate(psi_x_t**2, (x, -oo, oo)) - 1, A))
# {-sqrt(λ)*exp(I*t*ω), sqrt(λ)*exp(I*t*ω)}

Upvotes: 1

Related Questions