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

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[[1]];
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[[1]]["Grid"], grid}];
asol[t_, y_] = a[t, y] /. a -> bsol[[1]];]


Spatial domain and initial condition (IC):

Clear[g0, Ly, Uy]
Ly = -20; Uy = 20;
(*Giving a random IC for feq*)
SetAttributes[g0, Listable];
SeedRandom[1];
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[[1]] = 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[[1]] 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.

• ……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/… Commented Dec 20, 2018 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}]. Commented Dec 20, 2018 at 16:20
• @xzczd can 'pdetoode' handle two coupled PDEs? the PDEs above are decoupled. Commented Oct 30, 2023 at 21:31
• @S.Maths Of course. Notice pdetoode/pdetoae is just auxiliary tool for implementation of FDM, so it can of course handle PDE system, see e.g. mathematica.stackexchange.com/a/184285/1871 mathematica.stackexchange.com/q/278430/1871 mathematica.stackexchange.com/q/183745/1871 And you'll find more by searching pdetoode/pdetoae in this site. Commented Oct 31, 2023 at 1:41
• @xzczd thanks. the 1st link is close to my case. I adapted accordingly but still got the error "There are more dependent variables..." Commented Oct 31, 2023 at 10:10

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[[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,
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[
DeveloperToPackedArray@#["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[
DeveloperToPackedArray@#["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"]}


• Thank you. I just saw the answer. You didn't use pdetoae, so the related functions can be deleted. Commented Dec 21, 2018 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. Commented Dec 21, 2018 at 11:43