Skip to main content
1 of 3
Michael E2
  • 244.8k
  • 18
  • 351
  • 774

Here is another method that is more general in one way and less in another. It uses the values of the time derivative as the measure of when equilibrium has been reached. OTOH, it uses the MonitorMethod (see below), which does not work with multistep time-integration methods such the default LSODA.

The OP's example follows. It gives an error because, if there is a way with this method to stop the integration without causing NDSolve to give an error, it undocumented and I couldn't guess it. The feature of the MonitorMethod is that you can pass a "monitor function" to be called at each step (with access to any of the NDSolve state data). Here I tweaked it so that if the function returns False or $Failed, integration is stopped.

Clear[norm];
mf = Function[{h, sd, state, meth}, (* monitor function *)
   If[norm[ (* scaled norm of the time-derivative X' relative to X *)
      NDSolve`SolutionDataComponent[sd, "X'"],
      NDSolve`SolutionDataComponent[sd, "X"]
      ] < 1,
    Print["Reached equilibrium at t = ", 
     NDSolve`SolutionDataComponent[sd, "T"]];
    $Failed,
    True, True]];
usol = NDSolve[{pde, ic, bc}, uind, {t, 0, T}, {x, 0, 1}, 
  MaxStepSize -> 1,
  Method -> {MonitorMethod,
    "MonitorFunction" -> mf,
    "SaveNormAs" :> norm}]

MonitorMethod:

I added an option that allows the user to save the norm used by NDSolve in an external variable. It is used above to detect when the derivative is close enough to zero.

MonitorMethod // ClearAll;
MonitorMethod // Options = {Method -> Automatic, 
   "MonitorFunction" -> 
    Function[{h, sd, state, meth}, Print[{"H" -> h, "SD" -> sd}]],
   "SaveNormAs" -> Automatic};
MonitorMethod /: 
  NDSolve`InitializeMethod[MonitorMethod, stepmode_, sd_, rhs_, 
   state_, OptionsPattern[MonitorMethod]] := 
  Module[{submethod, mf, norm},
   mf = OptionValue["MonitorFunction"];
   submethod = OptionValue[Method];
   If[submethod === Automatic, submethod = "StiffnessSwitching"];
   submethod = 
    NDSolve`InitializeSubmethod[MonitorMethod, submethod, stepmode, 
     sd, rhs, state];
   norm = OptionValue[Automatic, Automatic, "SaveNormAs", Hold];
   norm /. {
     Hold[Automatic] :> {},
     Hold[nf_] :> (nf = state@"Norm")}; (* get the norm for the user *)
   MonitorMethod[submethod, mf]];
MonitorMethod[submethod_, mf_]["Step"[f_, h_, sd_, state_]] := 
  Module[{res},
   res = NDSolve`InvokeMethod[submethod, f, h, sd, state];
   If[Head[res] === NDSolve`InvokeMethod,
    Return[$Failed]];(* submethod not valid for monitoring *)
   If[! FreeQ[mf[h, sd, state, submethod], $Failed | False],
    (* mf indicates stopping *)
    Return[{0, $Failed, MonitorMethod[res[[-1]], mf]}; 
     "StopIntegration"](* causes error *)
    ];
   If[SameQ[res[[-1]], submethod],
    res[[-1]] = Null,
    res[[-1]] = MonitorMethod[res[[-1]], mf]];
   res];
MonitorMethod[___]["StepInput"] = {"Function"[All], "H", 
   "SolutionData", "StateData"};
MonitorMethod[___]["StepOutput"] = {"H", "SD", "MethodObject"};
MonitorMethod[submethod_, ___][prop_] := submethod[prop];

Other uses of MonitorMethod: (199228), (210480)

Michael E2
  • 244.8k
  • 18
  • 351
  • 774