0
$\begingroup$

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++;];
$\endgroup$
4
  • 2
    $\begingroup$ 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. ) $\endgroup$
    – xzczd
    Dec 20, 2018 at 12:53
  • $\begingroup$ ……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}] $\endgroup$
    – xzczd
    Dec 20, 2018 at 15:07
  • $\begingroup$ Dear @xzczd Thank you for the useful examples. $\endgroup$
    – user55777
    Dec 20, 2018 at 16:59
  • $\begingroup$ @user55777 If you explain correctly what you want to iterate, I will build a loop using debugged code. $\endgroup$ Dec 20, 2018 at 17:37

1 Answer 1

2
$\begingroup$

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

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.