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