2
$\begingroup$

Maybe it is a syntax problem, since I haven't been using Mathematica for much time, but I haven't been able to find a similar example in internet.

I am just trying to solve analitically the Sturm-Liouville problem of Schrödinger equation for a particle in a 2D box $[0,1 ] \times [0,1]$. Id est:

$ - \frac{\hbar^2}{2m} \Delta \psi(x,y) = E \psi(x, y) $

$\psi(0,y) = \psi(1,y) = \psi(x,0) = \psi(x,1) = 0$

This is what I wrote:

eqn = 2*pi^2*z[x, y] + (D[z[x, y], {x, 2}] + D[z[x, y], {y, 2}]) == 0;
bc = {z[x, 0] == 0, z[x, 1] == 0, z[0, y] == 0, z[1, y] == 0};
sol = DSolve[{eqn, bc}, z[x, y], {x, y}]

And this is what I get:enter image description here

$\endgroup$
  • 1
    $\begingroup$ You must have the equation and the initial/boundary conditions in a list. This however does not mean that DSolve can solve this equation. That would be a separate issue ! $\endgroup$ – Lotus Oct 5 '17 at 12:22
  • $\begingroup$ If you and @Miguel Bolín are the same person, please contact moderator for merging the accounts. $\endgroup$ – Alexey Popkov Oct 5 '17 at 12:43
2
$\begingroup$

As mentioned in the comments, the equation and boundary conditions need to be in the same list.

Beyond that, Mathematica has a little bit of difficulty solving the exact equation you've given: but if you're happy to sacrifice a bit of exactitude in the name of getting the right answer, you can try

eqn = (1 + $MachineEpsilon)*2*Pi^2*z[x, y] + (D[z[x, y], {x, 2}] + D[z[x, y], {y, 2}]) == 0; 
bc = {z[x, 0] == 0, z[x, 1] == 0, z[0, y] == 0, z[1, y] == 0}; 
sol = z[x, y] /. DSolve[Join[{eqn}, bc], z[x, y], {x, 0, 1}, {y, 0, 1}]

which gives

{2*Inactive[Sum][0, {K[1], 1, Infinity}]}

i.e., 0.

$\endgroup$
2
$\begingroup$

You can do this with M, but you need to help it. Start with the pde

 pde = D[\[Psi][x, y], x, x] + D[\[Psi][x, y], y, y] == 
-2 m En \[Psi][x, y]/\[HBar]^2

Separate variables

\[Psi][x_, y_] = X[x] Y[y]

pde = pde/\[Psi][x, y] // Expand

D[X[x], x, x]/X[x] + D[Y[y], y, y]/Y[y] == -((2 En m)/\[HBar]^2)

The first term is a fn of x only, the second a function of y and the rhs is a constant, so each term must be equal to a constant. A negative constant will give the desired sinusoidal solutions.

DSolve[{pde[[1, 1]] == -a^2, X[0] == 0}, X[x], x] // Flatten

{X[x] -> C[2] Sin[a x]}

Only put in one condition because, if not, it will return only the trivial solution of 0. To solve for the condition at x = 1 for a

X[x_] = X[x] /. % /. C[2] -> c1

Reduce[X[1] == 0, a]

    (C[1] \[Element] Integers && (a == 2 \[Pi] C[1] || a == 2 \[Pi] C[1] + \[Pi])) || c1 == 0

which boils down to

a = n Pi

and

$Assumptions = n \[Element] Integers

Now the y part. Subst the x part into the pde.

pde

D[Y[x], y, y]/Y[y] - n^2 Pi^2 == -((2 En m)/\[HBar]^2)

Set the Y term to another negative constant.

DSolve[{D[Y[y], y, y]/Y[y] == -b^2, Y[0] == 0}, Y[y], y] // Flatten

{Y[y] -> C[2] Sin[b y]}

Y[y_] = Y[y] /. % /. C[2] -> 1

Sin[b y]

I set the C[2] to 1 because it is multiplied by the constant in the x equation leaving only one constant anyway. The condition at y = 1

Y[1] == 0

Sin[b] == 0

and as before

b = k Pi

$Assumptions = $Assumptions && k \[Element] Integers

pde

-(\[Pi]^2 k^2) - \[Pi]^2 n^2 == -((2 En m)/\[HBar]^2)

Solve for the Energies

Solve[pde, En] // Flatten

{En -> (\[HBar]^2 (\[Pi]^2 k^2 + \[Pi]^2 n^2))/(2 m)}

Normalize the wave function

normint = Integrate[\[Psi][x, y]^2, {x, 0, 1}, {y, 0, 1}] == 1

c1^2/4 == 1

c1 = c1 /. Solve[normint, c1][[2]]

\[Psi][x, y]

2 Sin[\[Pi] k y] Sin[\[Pi] n x]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.