# A is fast, B is fast, but together they're Mathematica-crashing slow?

I'm trying to do something with finding solutions to a quantum mechanics problem with n wells. If there are 40 wells, I need to find the solution to an equation in the form:

{a,b}.MatrixPower[A,40].{c,d}


So, I use FindRoot and everything works out great, in a decent amount of time. The actual code used is something along these lines:

lowPropagator[
dist_, α_] :=  {{Cos[α dist],Sin[α dist]/α},
{-α Sin[α dist], Cos[α dist]}};
highPropagator[
dist_, β_] := {{Cosh[β dist], Sinh[β dist]/β},
{β Sinh[β dist],Cosh[β dist]}};
root[d_, n_, start_] :=
FindRoot[{α, 1}.
MatrixPower[lowPropagator[d, α].highPropagator[d,α], n - 1].
lowPropagator[d, \[Alpha]].{1,α},
{α,start}, MaxIterations -> 50
];
Manipulate[
root[d, n,α0], {{d, 1}, 0, 5}, {{n, 2}, 1, 10,
1}, {{α0, 1}, 0, 5}]


It's responsive even for larger n. Then, I try to plot the equation. The plotting function is a large/complicated piecewise function, alternating between the form $f=A \cos(\alpha x)+B\sin(\alpha x)$, and $f=A \cosh(\beta x)+B\sinh(\beta x)$ (where beta is a function of alpha). The matrices "lowPropagator" and "highPropagator" allow me to find the function and its derivative at each point where the form of the piecewise function changes, so I can use a matrix inverse to solve for A and B, and then plug all of this into a giant Piecewise[]. So there's a lot of computation going on in the function. The function is in the form:

getWaveFunction[α_, β_, \[Rho]_, a_, b_, n_] := Module[{...},
...
totalFunction[x_] :=
Piecewise[Table[{eqlist[[i]][x], x < wavePositions[[i]]}, {i, 1, 2 n}],
eqlist[[2 n + 1]][x]];
totalFunction[#] &
];


So the function isn't trivial, but it still runs fast, faster than the root finding even, especially when n is high. (it's doing at least 2n matrix multiplications)

However, when I do something like this:

plotIt[d_, n_, start_] := Module[{αnear, f},
αnear = root[d, n, start];
f = getWaveFunction[αnear, ..., d, n];
Plot[f[x], {x, -5, 5}]]


For some values, the function takes forever to draw the graph. The black evaluation bar doesn't pop up, so alt+. doesn't do anything. Basically, if I choose the input values incorrectly, most of the time Mathematica becomes completely unresponsive and I have to restart. I think Mathematica IS still running and calculating, because sometimes after 30 seconds or so it will return a result. Out of curiosity I left it running overnight with a different value and it was still frozen when I came back, so, something is not right.

So the biggest issue is that these two functions are fast on their own, at ALL input values, but at certain input values the combination of the two takes forever. My question is, how can I fix this to make the function run at a decent speed, and if I can't, how can I make it fail gracefully so that the notebook doesn't crash?

Relevant things:

• It fails at large n and large d or "dist".

• It IS taking a matrix with hyperbolic functions in it, to the nth power. So, the values can get pretty huge (understatement), and this is exactly what happens with the values at the crash. However, the root finding function on its own doesn't seem to crash at these values. Neither does the plotting function. It's only both together that causes problems.

• I figure evaluation might be the problem, maybe the root and/or the graph are trying to be found many times inside the plot function, but adding "//Evaluate" to the end of whatever f and αnear equal doesn't solve the problem.

Let me know if I should post the full code. The relevant bits add up to ~180 lines, with a good amount of spacing/comments. The full code produces graphics like this:

(The red line corresponds to the α result of FindRoot. The purpose of FindRoot is to find all solutions that are exponentially decaying on both sides of the graph, so the red line is exponentially decaying on both sides, but the dashed line [which corresponds to the "start" value in FindRoot] isn't a solution, since it diverges on the right side. I'm trying to find solutions to the time independent schrodinger equation, for a single particle, with n wells, and this is for the end of chapter 2 of "A First Course in Computational Physics" by DeVries [not homework, & the book doesn't ask for solutions this general].)

## Update with full code

So, I haven't been able to simplify the code to much significant extent. If I simplify it too much, to the point where it seems to work but isn't calculating what I want it to calculate, the problem seems to disappear! I do think the problem could be with evaluation though. If you comment out the "Show" at the end and just write "f2" in its place, (f2 is the function that is a result of getWaveFunction, with FindRoot passed in as an argument), it prints out, "totalFunction\$30552[#1] &" (or some similar number on the end). Since totalFunction doesn't exist in the scope it's at... I'm guessing that this is an evaluation problem, where Plot re-evaluates something it doesn't really have to a bunch of times.

Code:

(And, keep in mind while running it, it can crash)

(* \[Alpha] is the frequency of open regions
\[Beta] is the coefficient of high energy regions (should be \
\[Beta]=Sqrt[vm2h-\[Alpha]^2])
\[Rho] is the height at the leftmost wall
a is the width of the open region
b is the width of the high energy regions
n is the number of wells.

This function returns a pure function, that is supposed to be the \
wavefunction \[Psi][x]. (does not depend on time)
*)
getWaveFunction[\[Alpha]_, \[Beta]_, \[Rho]_, a_, b_, n_] :=
Module[{left, eqLeft, eqRight, eqlist, lowPropagator,
highPropagator, lowInverseCoefficients, highInverseCoefficients,
waveRecursiveFunction,
waveValues,
wavePositions,
waveConstants,
totalFunction},

lowPropagator[
dist_] :=  {{Cos[\[Alpha] dist],
Sin[\[Alpha] dist]/\[Alpha]},   {-\[Alpha] Sin[\[Alpha] dist],
Cos[\[Alpha] dist]}};
highPropagator[
dist_] := {{Cosh[\[Beta] dist],
Sinh[\[Beta] dist]/\[Beta]}, {\[Beta] Sinh[\[Beta] dist],
Cosh[\[Beta] dist]}};

(* The wavefunction and its derivative can be represented as:
A.{aconst,bconst}={\[Psi],\[Psi]'},
so to solve for aconst and bconst you just take the inverse of the \
matrix A,
where A either corresponds to the low solution or the high solution.
*)
lowInverseCoefficients[
x_] := {{Cos[\[Alpha] x], -Sin[\[Alpha] x]/\[Alpha]}, {Sin[\
\[Alpha] x], Cos[\[Alpha] x]/\[Alpha]}};
highInverseCoefficients[
x_] := {{Cosh[\[Beta] x], -Sinh[\[Beta] x]/\[Beta]}, {-Sinh[\
\[Beta] x], Cosh[\[Beta] x]/\[Beta]}};

(* Define the leftmost point so that everything is centered around \
0 *)
left = -(n*(a + b) - b)/2;

(* waveValues is to be an array giving the values of the \
wavefunction an its derivative. ie:
waveValues[[1]]=={\[Psi][left],\[Psi]'[left]};
waveValues[[2]]=={\[Psi][left+\[CapitalDelta]x],\[Psi]'[
left+\[CapitalDelta]x]};
etc.

wavePositions is just the position of the nth change in height,
with wavePositions[[1]] being the start of the leftmost dip.

waveConstants holds the values of coefficients to equations. So,
you'll have each waveConstants[[
n]] be equal to some vector {aconst, bconst}.
Whether the equation is in the zero energy region or high energy \
region depends on whether n is odd or even respectively.
So if n is odd the constants correspond to the equation (aconst*
Cos[\[Alpha] x]+bconst*Sin[\[Alpha] x]),
and if n is even the constants correspond to the equation (aconst*
Cosh[\[Beta] x]+bconst*Sinh[\[Beta] x]).
*)
waveValues = wavePositions = waveConstants = Table[0, {q, 1, 2 n}];
waveRecursiveFunction[1] := (
waveValues[[1]] = {\[Rho], \[Rho] \[Beta]};
wavePositions[[1]] = left;
waveConstants[[1]] =
lowInverseCoefficients[left].{\[Rho], \[Rho] \[Beta]};
);
waveRecursiveFunction[j_] := (
(* fill out previous values *)
waveRecursiveFunction[j - 1];
If[Mod[j, 2] == 0,
(* j=2m,
so to the left of it is a zero region and to the right is a \
high energy region. *)
waveValues[[j]] = lowPropagator[a].waveValues[[j - 1]];
wavePositions[[j]] = wavePositions[[j - 1]] + a; (*
width of the high energy region is a. *)
waveConstants[[j]] = highInverseCoefficients[0].waveValues[[j]];
,
(* j=2m+1,
so to the left of it is a high energy region and to the right \
is a zero energy region. *)
waveValues[[j]] = highPropagator[b].waveValues[[j - 1]];
wavePositions[[j]] = wavePositions[[j - 1]] + b;
waveConstants[[j]] =
lowInverseCoefficients[wavePositions[[j]]].waveValues[[j]];
];
);
waveRecursiveFunction[2 n];

eqlist =
Join[{Evaluate[(highInverseCoefficients[
0].{\[Rho], \[Rho] \[Beta]}).{Cosh[\[Beta] (# - left)],
Sinh[\[Beta] (# - left)]}] & },
Table[If[Mod[m, 2] == 0,
(Evaluate[
waveConstants[[m]].{Cosh[\[Beta] (# - wavePositions[[m]])],
Sinh[\[Beta] (# - wavePositions[[m]])]}]) &,
(Evaluate[
waveConstants[[m]].{Cos[\[Alpha] #], Sin[\[Alpha] #]}]) &
], {m, 1, 2 n}
]
];

totalFunction[x_] :=
Piecewise[
Table[{eqlist[[i]][x], x < wavePositions[[i]]}, {i, 1, 2 n}],
eqlist[[2 n + 1]][x]];
totalFunction[#] &
(*Show[Plot[totalFunction[x],{x,left-1,-left+1},PlotRange->{-1,1}],
drawPotential[a,b,1,n]]*)
(*{Plot[totalFunction[x],{x,-4,4}],waveValues,wavePositions,
waveConstants,eqlist,totalFunction[x]}//Column*)

];
Manipulate[
Module[{f1, f2, \[Alpha]near, left, lowPropagator, highPropagator},
c = Rationalize[c];
d = Rationalize[d];
a = Rationalize[a];
\[Rho] = Rationalize[\[Rho]];
left = -(n*(a + d) - d)/2;
f1 = getWaveFunction[\[Alpha], Sqrt[c - \[Alpha]^2], \[Rho], a, d,
n];
lowPropagator[
dist_, \[Alpha]_] :=  {{Cos[\[Alpha] dist],
Sin[\[Alpha] dist]/\[Alpha]},   {-\[Alpha] Sin[\[Alpha] dist],
Cos[\[Alpha] dist]}};
highPropagator[
dist_, \[Beta]_] := {{Cosh[\[Beta] dist],
Sinh[\[Beta] dist]/\[Beta]}, {\[Beta] Sinh[\[Beta] dist],
Cosh[\[Beta] dist]}};
\[Alpha]near = FindRoot[{Sqrt[c - \[Alpha]v^2], 1}.
MatrixPower[
lowPropagator[a, \[Alpha]v].highPropagator[d,
Sqrt[c - \[Alpha]v^2]], n - 1]
.lowPropagator[a, \[Alpha]v].{1, Sqrt[c - \[Alpha]v^2]},
{\[Alpha]v, \[Alpha]}, MaxIterations -> 50] // Evaluate;

\[Alpha]near = (\[Alpha]v /. \[Alpha]near);
If[Abs[Im[\[Alpha]near]] < 0.00001, \[Alpha]near =
Re[\[Alpha]near], \[Alpha]near = 0];
If[\[Alpha]near > 0,
f2 = Evaluate[
getWaveFunction[\[Alpha]near, Sqrt[c - \[Alpha]near^2], \[Rho],
a, d, n]], f2 = 0 &];

Show[Plot[f1[x],{x,left-1,-left+1},PlotStyle->Dashed],
Plot[f2[x],{x,left-1,-left+1},PlotStyle->Directive[Thick,Red],
Exclusions->None],
PlotRange->{-1.5,1.5}
]

], {{\[Alpha], 2}, 0.001, Sqrt[c]}, {{c, 10}, 0.001, 20}, {{a, 1},
0.001, 10}, {{d, 0.2}, 0.001, 10}, {{\[Rho], 0.2}, 0, 1}, {n, 2, 10,
1}]

-
It would be better to post a minimal example showing the problem. On the other hand, your problem is precisely that this happens with these functions, at this number of wells, so I appreciate that you may not be able to do that. Perhaps add the full code and an example showing the problem at the end? (I doubt I will have time to do this this week, but perhaps someone else will--and I'd certainly like to help). –  acl Mar 3 '13 at 21:10
As a side remark, perhaps an important lesson here is that this kind of approach to solving problems using computers doesn't scale beyond simple examples... –  acl Mar 3 '13 at 21:12
@acl Thanks. updated it w/ full code. Could you check my theory at the end of the post? Seems like it might be a red flag but I'm not sure what I can do about it. –  NeuroFuzzy Mar 4 '13 at 0:41
Does using Plot[Evaluate@f[x], {x, -5, 5}] help? –  m_goldberg Mar 4 '13 at 3:26
@m_goldberg Nope :/ (that's, Evaluate@f2[x], in the Show command at the bottom) –  NeuroFuzzy Mar 4 '13 at 4:20
show 1 more comment