Skip to main content
Clarification
Source Link
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. OTOHI approximated the derivative in my other answer, but since the derivative is computed at each step by NDSolve anyway, it makes sense to try to use it. OTOH, it uses the MonitorMethod (see below), which does not work with multistep time-integration methods such the default LSODA.

Clear[norm];
mf = Function[{ (* monitor function *)
     h,         (*  next step size (UNUSED) *)
     sd,        (*  SolutionData structure *)
     state,     (*  StateData object *)
     meth},     (* monitor functionMethodObject data structure (UNUSED) *)
   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, (* not necessary *)
  Method -> {MonitorMethod,
    "MonitorFunction" -> mf,
    "SaveNormAs" :> norm}]

Note: The step size was restricted because the step size quickly increases to 10. Without the restriction, it overshoots the time at which we are nearly at equilibrium by a bit. If you just want to ensure that equilibrium has been reached, then there is no point in including it.

The MonitorMethod invokes a submethod, but only certain ones are allowed. The documentation does not specify the disallowed methods, but the default Automatic time-integration method, LDSODA, turns out to be one them. Consequently, Automatic defaults to something else here, which I added an option thatchose to be "StiffnessSwitching". The method may be set with the Method suboption of MonitorMethod.

Other options include "SaveNormAs", which allows the user to save the norm used by NDSolve in an external variable, and "MonitorFunction", which sets the function to be called at each step. It The norm is constructed by NDSolve from the PrecisionGoal and AccuracyGoal. The norm is scaled so that norm < 1 indicates an error that is negligible. That criterion is used in mf[] above to detect when the derivative is close enough to zero. The arguments of the monitor function are mf[step_size, SolutionData, StateData, MethodObject]. To change them, change the call to mf[] found in the MonitorMethod[..]["Step"[..]] function.

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]];(* submethod not valid for monitoring *)
   (* call monitor function; stop integration if indicated *)
   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];

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.

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}]

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];

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. I approximated the derivative in my other answer, but since the derivative is computed at each step by NDSolve anyway, it makes sense to try to use it. OTOH, it uses the MonitorMethod (see below), which does not work with multistep time-integration methods such the default LSODA.

Clear[norm];
mf = Function[{ (* monitor function *)
     h,         (*  next step size (UNUSED) *)
     sd,        (*  SolutionData structure *)
     state,     (*  StateData object *)
     meth},     (*  MethodObject data structure (UNUSED) *)
   If[norm[ (* scaled norm of 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, (* not necessary *)
  Method -> {MonitorMethod,
    "MonitorFunction" -> mf,
    "SaveNormAs" :> norm}]

Note: The step size was restricted because the step size quickly increases to 10. Without the restriction, it overshoots the time at which we are nearly at equilibrium by a bit. If you just want to ensure that equilibrium has been reached, then there is no point in including it.

The MonitorMethod invokes a submethod, but only certain ones are allowed. The documentation does not specify the disallowed methods, but the default Automatic time-integration method, LDSODA, turns out to be one them. Consequently, Automatic defaults to something else here, which I chose to be "StiffnessSwitching". The method may be set with the Method suboption of MonitorMethod.

Other options include "SaveNormAs", which allows the user to save the norm used by NDSolve in an external variable, and "MonitorFunction", which sets the function to be called at each step. The norm is constructed by NDSolve from the PrecisionGoal and AccuracyGoal. The norm is scaled so that norm < 1 indicates an error that is negligible. That criterion is used in mf[] above to detect when the derivative is close enough to zero. The arguments of the monitor function are mf[step_size, SolutionData, StateData, MethodObject]. To change them, change the call to mf[] found in the MonitorMethod[..]["Step"[..]] function.

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 *)
   (* call monitor function; stop integration if indicated *)
   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];
added 3 characters in body
Source Link
Michael E2
  • 244.8k
  • 18
  • 351
  • 774

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 is 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.

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.

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 is 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.

Source Link
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)