2
$\begingroup$

I wrote a function that numerically solves a given Schrodinger equation and returns replacement rules to obtain the time evolution operator. It looks and is used as follows:

solver[H_, time_] := 
    soln = Module[{d, init, eqs, vars, solargs, t, t0, tf},
    d = Dimensions[H][[1]];
    t0 = time[[2]];
    tf = time[[3]];
    t = time[[1]];
    u[t_] := Table[Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}];
    init = Thread[Flatten /@ (u[t0] == IdentityMatrix[d])];
    eqs = Thread[Flatten /@ (I*u'[t] == H.u[t])];
    vars = Flatten[Table[Subscript[u, i, j], {i, 1, d}, {j, 1, d}]];
    solargs = Join[eqs, init];
    Return[NDSolve[solargs, vars, time, InterpolationOrder -> All, Method -> "ExplicitRungeKutta", AccuracyGoal -> 10]]
  ];

(* use like this, e.g. for 2D problem*)
ham[t_]:={{0,t},{t,0}};
solEvol=solver[ham[t],{t,0,20}];
evolOp[t_]=u[t]/.solEvol[[1]];

I do now need to work with the result obtained from my solver inside e.g. NMaximize to determine several parameters. To simplify things I reduced my original code/question to this snippet:

testFun[mat_] := Abs[mat];
uTemp2D[t_] := Table[Subscript[u, i, j], {i, 1, 2}, {j, 1, 2}];
hamParam[t_] = {{a, t}, {t, a}};
optimum = NMaximize[
    {testFun[uTemp2D[20] /. solver[hamParam[t], {t, 0, 20}]], 
    -10 < a < 10}, 
    a, 
    Method -> {"NelderMead"}
  ];

As far as I understand this now evaluates to

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.
ReplaceAll::reps: (* here come all entries of eqs in solver function*) is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing."

because NDSolve in solver does get non-numerical values (here a) as input. So the order of evaluation is incorrect. I thought of using something like _?NumericQ for matrices (https://mathematica.stackexchange.com/a/19600) but that is either working nor necessary accoring to a comment by MichaelE2.

Question

So in brief, what is the proper way to define my solver function or maybe the testFun to get results for my variable optimum?


Original code snippets (to which MichaelE2's answer refers)

Originally I was intending to execute the following (hParam[t] same as above):

partialTrace[states_, mat_] :=
  Module[{l = Length[states], trace},
   trace = Sum[{states[[i]]}\[Conjugate].mat.states[[i]], {i, 1, l}];
   Return[trace];
 ];

fidelPhase[evol_, Uideal_, states_] := 
  Module[{result}, 
   result = (1/(Length[states])^2)*(Abs[
   partialTrace[states, Uideal\[ConjugateTranspose].evol]])^2; 
   Return[result];];

solver[H_, time_, index_] := 
  soln = Module[{d, init, eqs, vars, solargs, t, t0, tf,meth = index}, 
  d = Dimensions[H][[1]];
  t0 = time[[2]];
  tf = time[[3]];
  t = time[[1]];
  u[t_] := Table[Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}];
  init = Thread[Flatten /@ (u[t0] == IdentityMatrix[d])];
  eqs = Thread[Flatten /@ (I*u'[t] == H.u[t])];
  vars = Flatten[Table[Subscript[u, i, j], {i, 1, d}, {j, 1, d}]];
  solargs = Join[eqs, init];
  Return[
    Which[
      meth == 1, NDSolve[solargs, vars, time, InterpolationOrder -> All, AccuracyGoal -> 10],
      meth == 2, NDSolve[solargs, vars, time, InterpolationOrder -> All, Method -> {"FixedStep", Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 15}}, StartingStepSize -> 1/300, AccuracyGoal -> 10],
      meth == 3, NDSolve[solargs, vars, time, InterpolationOrder -> All, Method -> "ExplicitRungeKutta", AccuracyGoal -> 10]
     ]
   ];
 ];

optimizer[gateTime_, ham_, ideal_, vars_, range_, states_] :=
  Module[{params},
  params = Map[Flatten, Transpose[{vars, range}]];
  solsOpt = NMaximize[
    fidelPhase[uTemp3D[gateTime]/.solver[ham/.tg -> gateTime, {t,0,gateTime},3][[1]], ideal,states],
  params,
  Method -> {"NelderMead", "Tolerance" -> Sqrt[$MachineEpsilon]}
  ];
  Return[solsOpt];
 ];

uIdeal = {{0, 1, 0}, {1, 0, 0}, {0, 0, 1}};
test = optimizer[
    20, hParam[t], uIdeal, 
    {a, b}, {{-10, 10}, {-10, 10}}, {{1, 0, 0}, {0, 1, 0}}
  ];

which now works with MichaelE2's answer.

$\endgroup$
7
  • 1
    $\begingroup$ solver[H_?(MatrixQ[#,NumericQ]&), time_, index_] :=.. requires your Hamiltonian consist of explicit numbers (no variables such as t, such as h[t] has). This should be fine, solver[H_, time_, index_] :=.., as you have in the beginning. At least it yields solutions. Since uTemp3D is not a variable of your system and not defined elsewhere, I don't know what to do with it or the rest of your problem. $\endgroup$ – Michael E2 Jan 27 '15 at 12:45
  • $\begingroup$ @MichaelE2 Indeed I forgot to add the uTemp3D defintion. I updated my question and added a concrete example of what is not working together with the error message(s) given. Hopefully the problem is a bit more clarified now. $\endgroup$ – Lukas Jan 27 '15 at 13:02
  • $\begingroup$ To get more and better answers, try to simplify your real problem to a minimal example. $\endgroup$ – Dr. belisarius Jan 27 '15 at 16:12
  • $\begingroup$ @belisarius I updated my question and reduced it to a minimal example for my problem. Hopefully it is now easier to read and understand. $\endgroup$ – Lukas Jan 28 '15 at 6:40
  • $\begingroup$ @MichaelE2 I updated my question and reduced it to a minimal example for my problem. Hopefully it is now easier to read and understand. $\endgroup$ – Lukas Jan 28 '15 at 6:40
3
$\begingroup$

Primary fixes: Adding First to fidelPhase and adding t to Subscript[u, i, j][t]. Making an objective function obj that is not evaluated until a and b are numeric may or may not be important. I won't have time to check it out (I've lost track of the original optimizer). The following works and it's not crucial to be so precise in one's fixes.

Clear[optimizer];
optimizer[gateTime_, ham_, ideal_, vars_, range_, states_] := 
 Module[{params, obj},
  obj[v_ /; VectorQ[v, NumericQ]] := 
   First@fidelPhase[
     uTemp3D[gateTime] /. 
      solver[ham /. Thread[vars -> v] /. tg -> gateTime, {t, 0, 
         gateTime}, 3][[1]], ideal, states];
  params = Map[Flatten, Transpose[{vars, range}]];
  solsOpt = 
   NMaximize[obj[vars], params, 
    Method -> {"NelderMead", "Tolerance" -> Sqrt[$MachineEpsilon]}, 
    EvaluationMonitor :> Print];
  Return[solsOpt];]

uTemp3D[t_] = Table[Subscript[u, i, j][t], {i, 1, 3}, {j, 1, 3}];
$\endgroup$
4
  • $\begingroup$ This works brilliant. Thanks alot! One further question: is this the approach one would usually choose for dealing with problems of that kind or is it just like the best fit to my specific problem? $\endgroup$ – Lukas Jan 28 '15 at 11:55
  • $\begingroup$ @Lukas I would say that coding the objective function as a separate obj with NumericQ protection of the NMaximize variables is a general technique, in the case where the function is defined in terms of numeric functions such as NDSolve, NIntegrate, FindRoot,.... (When the function is an algebraic formula, it's usually better not to have NumericQ to allow symbolic analysis.) Whether there's a better overall approach to your particular problem is not something I had time to think about.... $\endgroup$ – Michael E2 Jan 28 '15 at 13:38
  • $\begingroup$ ...I suspect the two minor "fixes" were hard to detect and confusing. One might try to think if there's a simpler way to build the solution, so that such errors are discovered earlier or more easily. $\endgroup$ – Michael E2 Jan 28 '15 at 13:41
  • $\begingroup$ Thanks for the quick answer and the brief explanation. I do highly appreciate any help and information to become better in using Mathematica. $\endgroup$ – Lukas Jan 28 '15 at 16:41

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.