Wie ersetze ich den Ausdruck im Nest?

Dec 11 2020

Ich muss die folgende Funktion mit wiederkehrendem Ersatz berechnen

f[x0_, y0_]:= (s1 + s2)/t1 /. {s1 -> 
NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0., 
   x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.}, 
  y, {t, 0, t1}][t1], s2 -> NDSolveValue[{x''[t] + x[t] == 0, y''[t] + x[t]^2 y[t] == 0.,
    x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 1.}, 
  y, {t, 0, t1}][t1]} /. {t1 -> Take[Reap[
    NDSolve[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0., 
      WhenEvent[x'[t] > 0., {Sow[t], "StopIntegration"}]}, 
     x, {t, 0., 100.}, 
     MaxStepSize -> 0.001]], {2, -1}][[1]][[1]][[1]]}

Aber ich treffe einen Fehler, der besagt "NDSolveValue: Endpunkt t1 in {t, 0., t1} ist keine reelle Zahl".

Der Schlüsselausdruck ist

NDSolveValue[{x''[t] + x[t] == 0, y''[t] + x[t]^2 y[t] == 0., 
 x[0.] == 1., x'[0.] == 0., y[0.] == 1., y'[0.] == 0.}, 
y, {t, 0, Re[t1]}][Re[t1]] /. {t1 -> Take[Reap[
    NDSolve[{x''[t] + x[t] == 0, x[0.] == 1., x'[0.] == 0., 
      WhenEvent[x'[t] > 0., {Sow[t], "StopIntegration"}]}, 
     x, {t, 0., 100.}, 
     MaxStepSize -> 0.001]], {2, -1}][[1]][[1]][[1]]}

Welches kommt der Fehler. Wie kann ich dieses Problem lösen?

Antworten

3 MichaelE2 Dec 12 2020 at 08:26

Nehmen Sie die t1Ersetzung in der s1,s2Ersetzung vor und gruppieren Sie sie in Klammern:

ff[x0_, y0_] := s1 + s2 /. (
    {s1 :> 
       NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0., 
          x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.}, 
         y, {t, 0, t1}][t1], 
      s2 :> NDSolveValue[{x''[t] + x[t] == 0, 
          y''[t] + x[t]^2 y[t] == 0., x[0.] == x0, x'[0.] == 0., 
          y[0.] == y0, y'[0.] == 1.}, y, {t, 0, t1}][t1]} /. {t1 ->
       NDSolveValue[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0., 
         WhenEvent[x'[t] > 0., {"StopIntegration"}]}, 
        Indexed[x["Domain"], {1, -1}], {t, 0., 100.}, 
        MaxStepSize -> 0.001]}
    );

ff[-1, -3]

(*  -2.54137  *)

Alternative:

ff[x0_, y0_] := 
  Block[{t1 = 
     NDSolveValue[{x''[t] + x[t] == 0., x[0.] == x0, x'[0.] == 0., 
       WhenEvent[x'[t] > 0., {"StopIntegration"}]}, 
      Indexed[x["Domain"], {1, -1}], {t, 0., 100.}, 
      MaxStepSize -> 0.001]},
   s1 + s2 /.
    {s1 :> 
      NDSolveValue[{x''[t] + x[t] == 0., y''[t] + x[t]^2 y[t] == 0., 
         x[0.] == x0, x'[0.] == 0., y[0.] == y0, y'[0.] == 0.}, 
        y, {t, 0, t1}][t1], 
     s2 :> NDSolveValue[{x''[t] + x[t] == 0, 
         y''[t] + x[t]^2 y[t] == 0., x[0.] == x0, x'[0.] == 0., 
         y[0.] == y0, y'[0.] == 1.}, y, {t, 0, t1}][t1]}
   ];