5
$\begingroup$

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.

Confusing Printed Information of pdetoode

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
$\endgroup$
5
  • $\begingroup$ Is pdetoode definition in general clear for you? $\endgroup$ Jan 13 at 2:15
  • $\begingroup$ @AlexTrounev Thanks for your attention. I think I have a clear understanding until the coordtoindex and the above delayed rule, and can't go further. Actually I tried to understand pdetoode by it application to this 3D time dependent pde here mathematica.stackexchange.com/questions/160012/… $\endgroup$
    – xinxin guo
    Jan 13 at 2:24
  • $\begingroup$ Do you understand that pdetoode is just implementation the method of lines with using NDSolve`FiniteDifferenceDerivative? $\endgroup$ Jan 13 at 2:36
  • $\begingroup$ @AlexTrounev Yes, I think I understand mol and I have carefully read through the "AdvancedNumericalDifferentialEquationSolvingInMathematica.pdf" several times. I can solve some pdes with Mathematica using 'mol'. But the code is ugly, not as general as pdetoode, and this is another movtivation that I started to try to understand 'pdetoode' $\endgroup$
    – xinxin guo
    Jan 13 at 2:42
  • $\begingroup$ Ok! Then the author of pdetoode has a good answer for you. $\endgroup$ Jan 13 at 6:28
5
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thanks for your help. In your response, I think I am clear before the last line "We also need to extract x, y, etc.". It is the pat that I am so confused. I thought dx1, dx2 are spatial derivative orders right? For example in Derivative[1,2,3,1][u][x,y,z,t], thenpat = Repeated[_, {3}], so dx1 = 1,2,3 should not mathch the pat. Thus, how can the spatial derivative orders be extracted by this rule? Could you please explain in more detail about this pattern? $\endgroup$
    – xinxin guo
    Jan 13 at 3:58
  • 2
    $\begingroup$ @xinxin No, it matches. Notice it's an unnamed pattern _. $\endgroup$
    – xzczd
    Jan 13 at 4:04
  • $\begingroup$ Oh, yeah! Finally! Your explanation is crystal clear! I mixed up Repeated[x_,{3}] and x: Reapeated[_,{3}] and the word Repeat in this command gives me a strong impression that something is repeated it self. Also, could please explain a little about coordtoindex function. It seems maps min and max coordinates to 1 and -1 and this makes sense to me. But, if I understand coorectly,why the other part that coordtoinex maps All to All? How the two Alls works in this function? $\endgroup$
    – xinxin guo
    Jan 13 at 4:33
  • $\begingroup$ @xinxin x: Repeated[_, {3}] still matches, but if it's Repeated[x_, {3}], then it won't match. As to coordtoindex, see my update. $\endgroup$
    – xzczd
    Jan 13 at 4:51
  • 1
    $\begingroup$ Thanks for your valuable help. The picture of pdetoode is becoming more and more clear for me with your explains. It would be take much longer for me to figure out "The function coordtoindex detects whether the discretized equation is a b.c. or not.". Actually, I and my classmates majored in civil engineering all learned so much from you in baidu tieba (we believe xzcyr is xzczd), mmaqa and here. We are not in famous 985 or 211 colleges, but we are so lucky to have guys like you in China, you makes us see deeply the power and beautify of Mathematica. Respect! $\endgroup$
    – xinxin guo
    Jan 13 at 5:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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