I understand mol but meet difficulty in understanding how `pdetoode` atumatically generate pde-to-ode-rules by using this strange pattern and rule?

Question

I often solve pdes for my research, and years ago I found pdetoode in this forum is very handy. Although it is a small piece of code, it solves several interesting and challenging pdes (here, here and here for example) in this forum.

Thus, I decied to figure out its working principle, because I strongly believe that it is greatly benificial to both ode and Mathematica understandings.

I have been studying hard on the package pdetoode by xzczd for many times and many days each time I picks it up.

I still can't figure out the how the pde converted to odes.

In particularlly, the following command:

((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])

I guess from the context that pat here means whatever repeats itself several times exactly. However, for code here, after I print dx1, and x1. In some cases dx1 = 0 1, x1=x y , this confues me a lot.

I also attempt to use Trace on pdetoode, but its long long out confuses me more.

How can this be? Maybe I misunderstand something here? What does this piece of code try to do?

I tried my best but I still fails to understand it.

Can anyone explain something on the above code or pdetoode? Thanks!

In order to analyze how pdetoode works, below is the code that I used to Print some local variables. The confusing printed information mentioned above is shown in the following figure.

Remove["Global`*"] // Quiet;

Clear[fdd, pde2ode, tooderule, pdetoae, rebuild]

(*====== fdd =======*)

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


(*====== pde2ode =======*)

pde2ode[funcvalue_List, rest__] := 
  With[{s = Style[#, Purple, Bold] &}, 
   Print[s@"pde2ode-1 ********************************"];
   Print[s@"funcvalue = ", s@funcvalue]; Print[s@"rest = ", s@rest];
   pde2ode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], 
    rest]];

pde2ode[{func__}[var__], rest__] := 
  With[{s = Style[#, Red, Bold] &}, 
   Print[s@"pde2ode-2 ********************************"];
   Print[s@"func = ", s@func]; Print[s@"var = ", s@var]; 
   Print[s@"rest = ", s@rest];
   pde2ode[Alternatives[func][var], rest]];

pde2ode[front__ (* front : u[x,t],t *), grid_?VectorQ, o_Integer, 
   periodic_ : False] := 
  With[{s = Style[#, Blue, Bold] &}, 
   Print[s@"pde2ode-3 ********************************"];
   Print[s@"front = ", s@front]; Print[s@"grid = ", s@grid];
   pde2ode[front, {grid}, o, periodic]];

pde2ode[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, PossibleZeroQ[# - #2[[1]]]}, {-1, 
             PossibleZeroQ[# - #2[[-1]]]}}, All] &, {coord, 
          bound}]]},
     With[{s = Style[#, Darker[Green], Bold] &}, 
      Print[s@"pde2ode-core *****************************"];
      Print[s@"func = ", s@func]; Print[s@"time = ", s@time]; 
      Print[s@"grid = ", s@grid];
      Print[s@"pos = ", s@pos]; Print[s@"bound = ", s@bound]; 
      Print[s@"spacevar = ", s@spacevar]; 
      Print[s@"coordtoindex = ", s@coordtoindex]];
     tooderule@Flatten@{
        (*------- 
        rule_1 --------*)
        ((u : func) | 
            Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat, 
          t_, 
          x2___] :>
         (With[{}, Print["-------------"]; 
           Print["u = ", u]; Print["dx1 = ", dx1]; 
           Print["pat = ", pat]; Print["dt = ", dt]; 
           Print["dx2 = ", dx2]; Print["x1 = ", x1]; 
           Print["x2 = ", x2]];
          Sow@coordtoindex@{x1, x2}; 
          With[{}, Print["coord = ", {x1, x2}]; 
           Print["coordtoindex[coord] = ", coordtoindex@{x1, x2}]];
          
          fdd[{dx1, dx2}, {grid}, 
           Outer[Derivative[dt][u@##]@t &, grid], 
           "DifferenceOrder" -> o, PeriodicInterpolation -> periodic]),
        (*------- rule_2 --------*)
        
        inde : spacevar :> 
         With[{i = Position[spacevar, inde][[1, 1]]}, 
          Outer[Slot@i &, grid]]
        }]]];


(*====== tooderule =======*)

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]


(*====== 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]] &@
   pde2ode[func[var, t], t, rest]]


(*====== rebuild =======*)

rebuild[funcarray_, grid_?VectorQ, timeposition_ : 1] := 
 rebuild[funcarray, {grid}, timeposition]

rebuild[funcarray_, grid_, timeposition_?Negative] := 
 rebuild[funcarray, grid, Range[Length@grid + 1][[timeposition]]]

rebuild[funcarray_, grid_, timeposition_ : 1] /; 
  Dimensions@funcarray === Length /@ grid := 
 With[{depth = Length@grid}, 
  ListInterpolation[
     Transpose[
      Map[Developer`ToPackedArray@#["ValuesOnGrid"] &, #, {depth}], 
      Append[Delete[Range[depth + 1], timeposition], timeposition]], 
     Insert[grid, Flatten[#][[1]]["Coordinates"][[1]], 
      timeposition]] &@funcarray]
(*=================== example of pdetoode ================*)

L = 1; T = 1; x0 = -L/4; sigma = L/30;
domain = {-L/2, L/2};
{eq1, eq2} = {D[G[x, y, t], t] == -(D[G[x, y, t], x] + D[G[x, y, t], y]) - 
     I (f[x, y, t] + f2[x, y, t]), 
   D[f[x, t], t] == -D[f[x, t], x] + f[x, t] - I (G2[x, t])};

{ic1, ic2} = {G[x, y, 0] == 
    Exp[-((x - x0)/(Sqrt[2] sigma))^2 - ((y - x0)/(Sqrt[2] sigma))^2], 
   f[x, 0] == 0};

{bc1, bc2} = {G[x, y, t] == 0 /. Outer[{# -> #2} &, {x, y}, domain], 
   f[x, t] == 0 /. List /@ Thread[x -> domain]};

points = 5;
grid = Array[# &, points, domain];

difforder = 2;

ptoofunc1 = 
 pde2ode[{G[x, y, t], f[x, y, t], f2[x, y, t]}, t, {grid, grid}, difforder]

grid
ptoofunc2 = pde2ode[{G2[x, t], f[x, t]}, t, grid, difforder];

ptoofunc1 = 
 pde2ode[{G, f, f2}[x, y, t], t, {grid, grid}, difforder]

ptoofunc2 = pde2ode[{G2, f}[x, t], t, grid, difforder];


del = #[[2 ;; -2]] &;
rule1 = {f[x_, y_][t_] :> f[x][t], f2[x_, y_][t_] :> f[y][t]};
ode1 = del /@ del@ptoofunc1@eq1;
ode1 = del /@ del@ptoofunc1@eq1 /. rule1;

rule2 = G2[x_][t_] :> G[x, 0][t];
ode2 = del@ptoofunc2@eq2 /. rule2;

odeic1 = ptoofunc1@ic1;

odeic2 = ptoofunc2@ic2;


diff = With[{sf = 1}, D[#, t] + #] &;
odebc1 = Map[diff, MapAt[del /@ # &, ptoofunc1@bc1, {1}], {-2}];

odebc2 = Map[diff, ptoofunc2@bc2, {-2}];

sol = NDSolveValue[{ode1, ode2, odeic1, odeic2, odebc1, 
    odebc2}, {Outer[G, grid, grid], f /@ grid}, {t, 0, T}];

solG = rebuild[sol[[1]], {grid, grid}, 3];
solf = rebuild[sol[[2]], grid, 2];

(*Manipulate[Plot3D[solG[x,y,t]//Evaluate,{x,##},{y,##},PlotRange->{-0.1,1}],{\
t,0,T}]&@@domain*)

Manipulate[
   Plot[solf[x, t] // Abs // Evaluate, {x, ##}, PlotRange -> {0, 0.2}], {t, 0,
     T}] & @@ domain

Answer 1

To understand the source code, let's first re-visit the syntax of pdetoode.

(* The following isn't a working sample,
   it's just for illustration. *)
pdetoode[u[variables], t, grid, differenceorder, periodic]
pdetoode[{u, v, …}[variables], t, {grid1, grid2, …}, differenceorder, periodic]
pdetoode[{u[variables], v[variables],…},t, {grid1, grid2, …}, differenceorder, periodic]

As mentioned there:

The syntax of pdetoode is as follows: 1st argument is the function to be discretized (which can be a list i.e. pdetoode can handle PDE system), 2nd argument is the independent variable in the resulting ODE system (usually it's the variable playing the role of "time" in the underlying model), 3rd argument is the list of spatial grid, 4th argument is difference order, 5th argument is to determine whether periodic b.c. should be set or not. (5th argument is optional, the default setting is False. )

Well, sometimes I feel that this syntax is a design miss, perhaps I should have made the syntax more simlar with NDSolve. We know, method of lines is a method that discretizes PDE(s) to a system of ODEs i.e. we need to discretize all the independent variables except for t (to be preciser, the independent variable of the resulting ODEs). For example, if the original dependent variable is u[x, t, y], it'll eventually be transformed to something equivalent to

(* The following isn't a working sample,
   it's just for illustration. *)
Outer[{x, y}|->u[x, y][t], gridx, gridy]

If it's Derivative[1, 1, 1][u][x, t, y], it's transformed to something equivalent to

(* The following isn't a working sample,
   it's just for illustration. *)
NDSolve`FiniteDifferenceDerivative[{1, 1}, {gridx, gridy}, 
                                   Outer[{x, y}|->u[x, y]'[t], gridx, gridy], …]

Notice this is done with a single fdd in the source code.

To automate the process, it's clear we need to know the position of t. This is programmatically detected by

(* The following isn't a working sample,
   it's just for illustration. *)
pos = Position[{var}, time][[1, 1]]

We also need to extract x, y, etc., because we need to determine whether we are discretizing b.c. or not (pdetoode can directly handle b.c. by design! ), and this is done by

(* The following isn't a working sample,
   it's just for illustration. *)
…
pat = Repeated[_, {pos - 1}]
…
… [x1 : pat, t_, x2___] :> (Sow@coordtoindex@{x1, x2}; …

Notice here the pattern inside Repeated is an unnamed pattern _, so it's just for restricting the length of the sequence, the elements of the sequence don't need to be exactly the same. If you still feel confused, try the following and think about why it works:

{a, b, c} /. {Repeated[_, {2}], c} :> "This matches of course"
(* "This matches of course" *)

The function coordtoindex detects whether the discretized equation is a b.c. or not. If it's left boundary, it outputs -1; right boundary, 1; not a boundary, All. These are all indexes for Part ([[]]). The output is throwed out by Sow and catched by Reap in the function

tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@ 
  Reap[expr /. rule]

For example, if we're discretizing u[x, t, y] in the domain $[1,4]\times[2,3]$, then u[1, t, y] will finally evaluate to something equivalent to

(* The following isn't a working sample,
   it's just for illustration. *)
Outer[{x, y}|->u[x, y][t], gridx, gridy][[-1, All]]

Hope now it is a bit clearer.