Solving two coupled PDEs with the pdetoode function posted on this site [closed]

I want to solve two coupled PDEs with pdetoode function:

feq = {D[ϕ, t] + U*D[ϕ, y] - D[ϕ, {y, 2}] - (s0 - 1/20*y^2)*ϕ == 0}

and

beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0}

To this end, I encapsulated them in the two solvers, which are very similar (function pdetoode not included here, please find it from the above link):

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

beqSolver[s_, T_, aT_] :=
Module[{s0 = s, Tend = T, beqIC = aT},
With[{a = a[t, y], U = y},
beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0};
bic = {a == beqIC[y]} /. t -> Tend;
bbc = {{a == 0} /. y -> Ly, {a == 0} /. y -> Uy};];
bptoofunc = pdetoode[a[t, y], t, grid, difforder];
bdel = #[[2 ;; -2]] &;
bode = bdel@bptoofunc@beq[];
bodeic = bptoofunc@bic;
bodebc = bptoofunc@bbc;
bsollst =
NDSolveValue[{bodebc, bodeic, bode}, Map[a, grid], {t, Tend, 0}];
bsol = ListInterpolation[
DeveloperToPackedArray@#["ValuesOnGrid"] & /@ bsollst //
Transpose, {Flatten@bsollst[]["Grid"], grid}];
asol[t_, y_] = a[t, y] /. a -> bsol[];]

Spatial domain and initial condition (IC):

Clear[g0, Ly, Uy]
Ly = -20; Uy = 20;
(*Giving a random IC for feq*)
SetAttributes[g0, Listable];
SeedRandom;
g0[y_?NumericQ] =
BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 39], {0.}],
SplineClosed -> False][(y + 20)/40];

Define two functions: J is a function of interest, aT is a function used for constructing the IC for beq

Clear[J, aT];
J[g_, phiT_, y_] :=
NIntegrate[Dot[g[y], g[y]], {y, Ly, Uy}]/
NIntegrate[Dot[phiT[y], phiT[y]], {y, Ly, Uy}];
(*define how to generate IC for beq from the soln of feq*)
aT[g_, phiT_,
y_] := (-2*phiT[y]*
NIntegrate[Dot[g[y], g[y]], {y, Ly, Uy}])/(NIntegrate[
Dot[phiT[y], phiT[y]], {y, Ly, Uy}])^2;

Solve the two equation:

s0 = 0.4; Tend = 10;
g[y_?NumericQ] = g0[y];
phisol[t, y] = feqSolver[s0, Tend, g];
phiT[y] = phisol[Tend, y];
Jtest[] = J[g, phiT, y]; (*evaluate a function of interest and save it as the first element of a list Jtest*)
aT[y] = aT[g, phiT, y];
asol[t, y] = beqSolver[s0, Tend, aT];
a0[y] = asol[0, y];

When run this code, MMA mainly gives this kind of error:

Part::take: "Cannot take positions 2 through -2 in pdetoode[ϕ[t,y],t,grid,difforder]" \$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

By running the code line by line, I found that the first error results from phisol[t, y] = feqSolver[s0, Tend, g];

In the above code, I have used Jtest[] since in my real problem I want to evaluate and save Jtest[[iter]] continuously (in a list) in a loop with an iterator iter. Then Jtest[[iter]] will be used for stopping the loop as follows:

Jtest[[iter]] = J[g, phiT, y]; AppendTo[Jtest, Jtest[[iter]]];
If[Abs[Jtest[[iter]] - Jtest[[iter - 1]]] < acc, Break[]];

**My problems are (1) should I define an empty list in the beginning, i.e. something like Jtest={}?; and (2) how do I fix these problems?

I have searched on the SE site for a similar problem but failed to solve it so far. And I have tried my best to learn the fundamental syntax of Mathematica for several months but I'm still in trouble. Please help and any suggestions and comments are very welcome. Thank you in advance.

closed as off-topic by xzczd, bobthechemist, bbgodfrey, mikado, SzabolcsJan 2 at 13:35

This question appears to be off-topic. The users who voted to close gave these specific reasons:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – bbgodfrey, mikado, Szabolcs
• "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – xzczd, bobthechemist
If this question can be reworded to fit the rules in the help center, please edit the question.

• ……You still don't understand why J[[iter]] = J[g, phiT, y]; doesn't make sense at all. Please consider asking a question about this, rather than immediately trying to deal with such a big problem. Also, recently I've written a relevant Chinese tutorial, please consider reading it: note.youdao.com/share/… – xzczd Dec 20 '18 at 10:21
• @user55777 Are you sure about this equation beq = {D[a, t] + U*D[a, y] + D[a, {y, 2}] + (s0 - 1/20*y^2)*a == 0} ? I ask about the sign + D[a, {y, 2}]. – Alex Trounev Dec 20 '18 at 16:20

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[
Developer`ToPackedArray@#["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"]} • Thank you. I just saw the answer. You didn't use pdetoae, so the related functions can be deleted. – user55777 Dec 21 '18 at 2:33
• @user55777 I just used your code in which I added parameters grid, difforder and correctly used pdetoode. There is no reference to 'pdetoae' in your code. – Alex Trounev Dec 21 '18 at 11:43