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