# How to use a test function correctly to terminate a While loop?

I define a test function in order to stop my While loop:

J[g_, phiT_, y_] := NIntegrate[Dot[g[y], g[y]], {y, Ly, Uy}]/NIntegrate[Dot[phiT[y], phiT[y]], {y, Ly, Uy}];


where Ly and Uy are constants to specify the spatial domain, both g[y] and phiT[y] are functions that are obtained within a loop.

In my code, I want to evaluate and save J[[iter]] continuously in the loop with an iterator iter. Then J[[iter]] will be used for stopping the While loop as follows:

While[True,
(*some code here*)
J[[iter]] = J[g, phiT, y]; AppendTo[J, J[[iter]]];
If[Abs[J[[iter]] - J[[iter - 1]]] < acc, Break[]];
(*some code here*)
iter++;];


As indicated by @xzczd, J[[iter]] = J[g, phiT, y] is meaningless. My understanding is as follows:

J[[iter]] means that one is going to access an element (indexed by iter) of a list, named J. However, J is the name of a fucntion defined in the beginning. This is an apparent contradiction.

But I cannot figure out how to resolve this problem. The main reason, I think, is that I do not know how many loops are needed and that's why I used the While loop.

In contrast, if I know how many loops are involved, then I can come up with, for example, the following code using Do loop:

ti = 0; tf = 10; \[Delta]t = 0.01;
ItrNo = Round[(tf - ti)/\[Delta]t];
Jtest = ConstantArray[0, ItrNo - 1];
Do[(*some code here*);
Jtest[[iter]] = J[g, phiT, y];
(*some code here*),
{iter, ItrNo - 1}]


Can anybody help me? Thank you in advance!

Here is a small example, though may not be minimal, since I want to use pdetoode. Please copy the function pdetoode from the above link.

feqSolver[s_, T_, g_] := Module[{s0 = s, Tend = T, feqIC = g},
With[{\[Phi] = \[Phi][t, y]},
feq = {D[\[Phi], t] + y*D[\[Phi], y] -
D[\[Phi], {y, 2}] - (s0 - 1/20*y^2)*\[Phi] == 0};
fic = {\[Phi] == feqIC[y]} /. t -> 0;
fbc = {{\[Phi] == 0} /. y -> Ly, {\[Phi] == 0} /. y -> Uy};];
fptoofunc = pdetoode[\[Phi][t, y], t, grid, difforder];
fdel = #[[2 ;; -2]] &;
fode = fdel@fptoofunc@feq[];
fodeic = fptoofunc@fic;
fodebc = fptoofunc@fbc;
fsollst = NDSolveValue[{fodebc, fodeic, fode},
Map[\[Phi], grid], {t, 0, Tend}];
fsol = ListInterpolation[DeveloperToPackedArray@#["ValuesOnGrid"] & /@ fsollst//Transpose, {Flatten@fsollst[]["Grid"], grid}];
phisol[t_, y_] = \[Phi][t, y] /. \[Phi] -> fsol[];]

s0 = 0.4; Tend = 5; acc = 10^-4;
Ly = -20; Uy = 20; iter = 1;
(*Random IC*)
SetAttributes[g0, Listable];
SeedRandom;
g0[y_?NumericQ] = BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 39], {0.}],
SplineClosed -> False][(y + 20)/40];
g[y] = g0[y];

J[g_, phiT_, y_] := g[y]/phiT[y];

While[True, gOld[y] = g[y];
phisol[t, y] = feqSolver[s0, Tend, gOld];
phiT[y] = phisol[Tend, y];
Jtest[[iter]] = J[g, phiT, y];
AppendTo[Jtest, Jtest[[iter]]];
If[Abs[Jtest[[iter]] - Jtest[[iter - 1]]] < acc, Break[]];
g[y] = -1/2*phiT[y];
iter++;];

• You should put more effort in creating a minimal example. The definition of J doesn't need to be that complicated. You can even simplify it to e.g. J[x_]=x^2. (The stopping criterion should be adjusted accordingly, of course. ) Dec 20, 2018 at 12:53
• ……Think carefully about the following questions. 1. How does this sample work?: Clear[J]; J[x_, y_] = x^2 + y; iter = 1; lst = {}; AppendTo[lst, J[iter, RandomReal[]]]; While[True, iter++; AppendTo[lst, J[iter, RandomReal[]]]; If[lst[[iter]] - lst[[iter - 1]] > 20, Break[]]]; lst; 2. How does this sample work?: Clear[J]; J[x_, y_] = x^2 + y; iter = 1; J[iter] = J[iter, RandomReal[]]; While[True, iter++; J[iter] = J[iter, RandomReal[]]; If[J[iter] - J[iter - 1] > 20, Break[]]]; Print@DownValues[J]; Table[J[i], {i, 1, iter}] Dec 20, 2018 at 15:07
• Dear @xzczd Thank you for the useful examples. Dec 20, 2018 at 16:59
• @user55777 If you explain correctly what you want to iterate, I will build a loop using debugged code. Dec 20, 2018 at 17:37

I debugged the code for the first iteration.

Clear[fdd, pdetoode, tooderule]
fdd[{}, grid_, value_, order_, periodic_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[front__, grid_?VectorQ, o_Integer, periodic_: False] :=
pdetoode[front, {grid}, o, periodic];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer,
periodic : True | False | {(True | False) ..} : False] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, # === #2[]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o,
PeriodicInterpolation -> periodic]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]
Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[pde, #[
pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]]

g = Interpolation[
Table[{-20 + i, If[1 < i < 40, RandomReal[{-1, 1}], 0]}, {i, 0,
40}], InterpolationOrder -> 2];
s0 = 0.4; Tend = 10;

feq = {D[f[y, t], t] + y*D[f[y, t], y] -
D[f[y, t], {y, 2}] - (s0 - 1/20*y^2)*f[y, t] == 0};
fic = {f[y, 0] == g[y]};
fbc = {f[Ly, t] == 0, f[Uy, t] == 0};
points = 25;
difforder = 2;
{Ly, Uy} = domain = {-20, 20};
grid = Array[# &, points, domain];
fptoofunc = pdetoode[f[y, t], t, grid, difforder];
fdel = #[[2 ;; -2]] &;
fode = fdel /@ fptoofunc@feq;
fodeic = fdel /@ fptoofunc@fic;
fodebc = fbc // fptoofunc;
fsollst =
NDSolveValue[{fodebc, fodeic, fode},
f /@ grid // Flatten, {t, 0, Tend}];
fsol = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ fsollst //
Transpose, {Flatten@fsollst[]["Grid"], grid}];

int1 = NIntegrate[g[y]^2, {y, Ly, Uy}];
int2 = NIntegrate[fsol[Tend, y], {y, Ly, Uy}];
h[y_] := -2*fsol[Tend, y]*int1/int2^2;
Jtest =
int1/int2;(*evaluate a function of interest and save it as the first \
element of a list Jtest*)
Clear[f];
feq = {D[f[y, t], t] + y*D[f[y, t], y] +
D[f[y, t], {y, 2}] - (s0 - 1/20*y^2)*f[y, t] == 0};
fic = {f[y, 0] == h[y]};
fbc = {f[Ly, t] == 0, f[Uy, t] == 0};
points = 25;
difforder = 2;
{Ly, Uy} = domain = {-20, 20};
grid = Array[# &, points, domain];
fptoofunc = pdetoode[f[y, t], t, grid, difforder];
fdel = #[[2 ;; -2]] &;
fode = fdel /@ fptoofunc@feq;
fodeic = fdel /@ fptoofunc@fic;
fodebc = fbc // fptoofunc;
fsollst =
NDSolveValue[{fodebc, fodeic, fode},
f /@ grid // Flatten, {t, 0, Tend}];
fsol1 = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ fsollst //
Transpose, {Flatten@fsollst[]["Grid"], grid}];

{Plot[g[y], {y, -20, 20}, AxesLabel -> {"y", "g"}],
Plot[h[y], {y, -20, 20}, AxesLabel -> {"y", "h"}],
ContourPlot[fsol[t, y], {t, 0, Tend}, {y, -20, 20}, Contours -> 20,
PlotRange -> All, PlotLegends -> Automatic, ColorFunction -> Hue,
FrameLabel -> {"t", "y"}, PlotLabel -> "\[Phi]"],
ContourPlot[fsol1[t, y], {t, 0, Tend}, {y, -20, 20}, Contours -> 20,
PlotRange -> All, PlotLegends -> Automatic, ColorFunction -> Hue,
FrameLabel -> {"t", "y"}, PlotLabel -> "a"]} 