0

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[[1]];
  fodeic = fptoofunc@fic;
  fodebc = fptoofunc@fbc;
  fsollst = NDSolveValue[{fodebc, fodeic, fode}, 
  Map[\[Phi], grid], {t, 0, Tend}];
  fsol = ListInterpolation[Developer`ToPackedArray@#["ValuesOnGrid"] & /@ fsollst//Transpose, {Flatten@fsollst[[1]]["Grid"], grid}];
  phisol[t_, y_] = \[Phi][t, y] /. \[Phi] -> fsol[[1]];]

s0 = 0.4; Tend = 5; acc = 10^-4;
Ly = -20; Uy = 20; iter = 1;
(*Random IC*)
SetAttributes[g0, Listable];
SeedRandom[1];
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++;];
user55777
  • 671
  • 3
  • 16
  • 2
    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. ) – xzczd Dec 20 '18 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}] – xzczd Dec 20 '18 at 15:07
  • Dear @xzczd Thank you for the useful examples. – user55777 Dec 20 '18 at 16:59
  • @user55777 If you explain correctly what you want to iterate, I will build a loop using debugged code. – Alex Trounev Dec 20 '18 at 17:37

1 Answers1

2

I debugged the code for the first iteration.

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

pdetoode[funcvalue_List, rest__] := 
  pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], 
   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, 
        MapThread[
         Piecewise[{{1, # === #2[[1]]}, {-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[[1]], 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[
   Developer`ToPackedArray@#["ValuesOnGrid"] & /@ fsollst // 
    Transpose, {Flatten@fsollst[[1]]["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[1] = 
 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[
   Developer`ToPackedArray@#["ValuesOnGrid"] & /@ fsollst // 
    Transpose, {Flatten@fsollst[[1]]["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"]}

fig1

Alex Trounev
  • 44,369
  • 3
  • 48
  • 106