Trouble to assign solution of differential equation to a function

Question

I have written a function to Solve the following differential equation enter image description here where the number n (called in my code Nphoton) is a variable. Nevertheless I have troubles to assign the solutions sol of the differential equation to functions funCa0, funCb0 ... The code

Do[
    ToExpression[
      "fun" <> ToString[vec[[i]]] <> "=Evaluate[" <> 
       ToString[vec[[i]]] <> "[t]/.sol][[1]]"];
    , {i, 1, Length[vec]}
    ];

returns the error message

ReplaceAll::reps: {sol} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

I do not understand this message, what is wrong with sol?

My entire function is (it can be run simply by pasting it):

funoscillation[Nphoton_, intialAtomState_, NphotonInitial_, g0_, g1_,wSP_, w0_] := 
  Module[{g0li, g1li, Delta, Omega, vec, funvec, eqns, sol, funPlot, 
    solReal, end},
   g0li = Conjugate[g0];
   g1li = Conjugate[g1];

   Delta = wSP - w0;
   Omega = wSP + w0;

   If[Nphoton - 1 < NphotonInitial, 
    Print["ATTENTION: no intial state provided"]; Goto[end];];
   vec = {};
   funvec = {};
Do[
  AppendTo[vec, ToExpression[StringJoin["Ca", ToString[i]]]];
  AppendTo[funvec, 
   ToExpression[StringJoin["funCa", ToString[i]]]]
  , {i, 0,    Nphoton - 1}];
Do[
  AppendTo[vec, ToExpression[StringJoin["Cb", ToString[i]]]];
  AppendTo[funvec, 
   ToExpression[StringJoin["funCb", ToString[i]]]]
  , {i, 0,    Nphoton - 1}];

   eqns = {};
   Do[
    If[i == 0,
     AppendTo[eqns, 
      ToExpression[
       "Ca" <> ToString[i] <> "'[t]==-I " <> ToString[g0] <> 
        " Sqrt[" <> ToString[i] <> 
        "+1] Exp[-I " <> ToString[Delta] <> " t] Cb" <> 
        ToString[i + 1] <> "[t]"]];

        AppendTo[eqns, 
          ToExpression[

       "Cb" <> ToString[i] <> "'[t]==-I " <> ToString[g1] <> 
        " Sqrt[" <> ToString[i] <> 
              "+1] Exp[-I " <> ToString[Omega] <> " t] Ca" <> 
        ToString[i + 1] <> "[t]"]];
        ];

      If[i > 0 && i < Nphoton - 1,
        AppendTo[eqns, 
          ToExpression[

       "Ca" <> ToString[i] <> "'[t]==-I " <> ToString[g0] <> 
        " Sqrt[" <> ToString[i] <> 
              "+1] Exp[-I " <> ToString[Delta] <> " t] Cb" <> 
        ToString[i + 1] <> 
              "[t]-I " <> ToString[g1li] <> " Sqrt[" <> ToString[i] <> 
        "] Exp[I " <> ToString[Omega] <> " t] Cb" <> 
              ToString[i - 1] <> "[t]"]];

        AppendTo[eqns, 
          ToExpression[

       "Cb" <> ToString[i] <> "'[t]==-I " <> ToString[g1] <> 
        " Sqrt[" <> ToString[i] <> 
              "+1] Exp[-I " <> ToString[Omega] <> " t] Ca" <> 
        ToString[i + 1] <> 
              "[t]-I " <> ToString[g0li] <> " Sqrt[" <> ToString[i] <> 
        "] Exp[I " <> ToString[Delta] <> " t] Ca" <> 
              ToString[i - 1] <> "[t]"]];
        ];

      If[i == Nphoton - 1,
        AppendTo[eqns, 
          ToExpression[

       "Ca" <> ToString[i] <> "'[t]==-I " <> ToString[g1li] <> 
        " Sqrt[" <> ToString[i] <> 
              "] Exp[I " <> ToString[Omega] <> " t] Cb" <> 
        ToString[i - 1] <> "[t]"]];

        AppendTo[eqns, 
          ToExpression[

       "Cb" <> ToString[i] <> "'[t]==-I " <> ToString[g0li] <> 
        " Sqrt[" <> ToString[i] <> 
              "] Exp[I " <> ToString[Delta] <> " t] Ca" <> 
        ToString[i - 1] <> "[t]"]];
        ];

      , {i, 0, Nphoton - 1}];


    Do[

      If[i == NphotonInitial ,
      If[intialAtomState == b,
       AppendTo[eqns, ToExpression["Ca" <> ToString[i] <> "[0]==0"]];
       AppendTo[eqns, ToExpression["Cb" <> ToString[i] <> "[0]==1"]];];
      If[intialAtomState == a,
       AppendTo[eqns, 
        ToExpression["Ca" <> ToString[i] <> "[0]==1"]];    
       AppendTo[eqns, ToExpression["Cb" <> ToString[i] <> "[0]==0"]];],

      AppendTo[eqns, ToExpression["Ca" <> ToString[i] <> "[0]==0"]];
      AppendTo[eqns, ToExpression["Cb" <> ToString[i] <> "[0]==0"]];
          ];

      , {i, 0, Nphoton - 1}];

   sol = NDSolve[eqns, vec, {t, 0, 10}];

   Do[
    ToExpression[
      "fun" <> ToString[vec[[i]]] <> "=Evaluate[" <> 
       ToString[vec[[i]]] <> "[t]/.sol][[1]]"];
    , {i, 1, Length[vec]}
    ];

   Print[funvec];
   Print[funCa0];
   funPlot = Abs[(# & /@ funvec)]^2;

   Plot[funPlot, {t, 0, 10}, PlotLegends -> vec, PlotRange -> All]

   ];

funoscillation[2, a, 1, 1*0, 1, 2 Pi, 2 Pi]

That's a lot of code to read for what is probably a rather simple question. Could you try to provide a more minimal example? Also, have you read this?: (6669) — Mr.Wizard, Sep 12 '14 at 11:32
Also, please see this question as well as the five links in the comment below it: (6511) -- you should not have to use ToString and ToExpression very often; here they are clearly being "abused" for a purpose that would be better served with a different method. — Mr.Wizard, Sep 12 '14 at 11:36
This is indeed a long code but the first part is only to define the differential equation eqns to solve. You can start reading almost at the end to understand my problem: from sol = NDSolve[eqns, vec, {t, 0, 10}]; (the last ten lines) — lambertmular, Sep 12 '14 at 12:02
So why not use a shorter example, e.g. sol = NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t]^3, x[0] == y[0] == 1}, {x, y}, {t, 20}] (from the documentation) unless it lacks a certain property of your actual code? — Mr.Wizard, Sep 12 '14 at 12:04
Yes it is different. I edited my post. The problem comes from the fact that the number of differential equations to solve is not fixed. — lambertmular, Sep 12 '14 at 12:17
I get the error "ReplaceAll::reps: {sol} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>", which would be helpful to include in the question. Such an error ought to be investigated, and it is a strong hint to experienced users what to look for. Not everyone is going to paste your code into M and see what happens. Having a hint like that in the question might encourage others to investigate, which is what you want to do. Also, if you were to put Print statements between the Do loops, you could localize the problem. — Michael E2, Sep 12 '14 at 14:07
Sorry, I forgot to mention it. I erase the Print statements before posting, next time I will conserve them. — lambertmular, Sep 12 '14 at 15:26

Michael E2 · Answer 1 · 2014-09-12T13:30:44.473

I realize from the idiomatic way this program was written that Mathematica is not a language you're very familiar with yet, at least not with the efficient and readable ways Mathematica can address your project. I won't try to suggest every possible improvement but focus on the specific problem of getting the functions you want into the Plot.

The following does not work because the string "sol" is converted to Global`sol and does not match the Module variable sol, which is localized in the form sol$nnn where nnn is some integer you don't really have control over. The other problem is that funvec is not updated here, which is what I assumed you wanted to have happen. (funvec has the value {} when you try to construct funPlot at the end.) Instead, different variables are created and never used.

Do[ToExpression[
    "fun" <> ToString[vec[[i]]] <> "=Evaluate[" <> 
     ToString[vec[[i]]] <> "[t]/.sol][[1]]"];, {i, 1, Length[vec]}];

I would advise a different approach. (Edit: fixed typos/misunderstandings of OP's code.)

funvec = Flatten@
     Table[ToExpression[pre <> ToString[i]], {i, 0, Nphoton - 1}, {pre, {"Ca", "Cb"}}];

funPlot = Abs[Through[funvec[t]] /. sol]^2;

Plot[funPlot, {t, 0, 10}, PlotLegends -> funvec, PlotRange -> All, 
   Axes -> False, Frame -> True]];

Mathematica graphics

(Ca0[t] and Cb1[t] seem to be identically zero.)

In the above, funvec is the list of symbols {Ca0, Cb0, Ca1, Cb1}. You could probably use this elsewhere in the code.

Theoretically one could use Block[{sol}, ...] instead of Module for the variable sol, and the Do loop would evaluate properly, if it created the list funvec. Note that your Do loops are used to create lists. This is usually done with Table in the form I used above for funvec. (Flatten is not always necessary.)

I want to assign a function to every solution of my differential equation (Ca0,Cb0,...). This is what I wanted to do with the code Do[ ToExpression[ "fun" <> ToString[vec[[i]]] <> "=Evaluate[" <> ToString[vec[[i]]] <> "[t]/.sol][[1]]"]; , {i, 1, Length[vec]} ];. Your code is not doing it — lambertmular, Sep 12 '14 at 16:02
I have edited my code so now funvec is well defined. But the problem comes from sol — lambertmular, Sep 12 '14 at 16:13
@lambertmular Quite so: The problem with sol is described in my answer. Re first comment: Your code doesn't work, so it was hard to figure out what to do with it. It's a very idiomatic way to program Mathematica. You can get the InterpolatingFunctions with {funCa0, funCb0,...} = funvec /. First[sol], constructing the LHS as you will. Your main problem is using sol in a string with ToExpression, when sol is Module variable. — Michael E2, Sep 12 '14 at 18:47

Trouble to assign solution of differential equation to a function

1 Answers1