# How to resolve the fail when applying a user-defined solver?

The user-defined function pdetoode and pdetoae, developed by @xzczd, is very useful to deal with a PDE system when there is difficult to discretize it into a system of ODEs. In this situation, NDSolve often transforms it into a DAE system and then uses a DAE solver, which is weaker than ODE solver. That is why it normally failed in such a case.

Here I tried to use pdetoae to solve such a PDE system with random initial conditions. Then I want to use the solutions to solve for another equation further by calling the previous solutions somehow. Please find the function pdetoode and pdetoae at this link.

Clear[f, m, Tend];
f[y_] := y; m = 300; Tend = 10;

(*Ramdom Iintial Conditions*)
Clear[a0, b0, c0, d0]
SeedRandom;
a0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
b0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
c0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
d0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];

Clear[difforder, domain, points, grid];

difforder = 4;
domain[y] = {-1, 1}; domain[t] = {0, Tend};
points[y] = 201; points[t] = 101;
(grid@# = Array[# &, points@#, domain@#]) & /@ {y, t};

Clear[a, b, c, d];
With[{a = a[t, y], b = b[t, y], c = c[t, y], d = d[t, y]},
eq = {I*a + D[b, y] + I*c == 0,
D[a, t] + I*f[y]*a + b*D[f[y], y] + I*d ==
1/m*(D[a, {y, 2}] - 2*a),
D[b, t] + I*f[y]*b + D[d, y] == 1/m*(D[b, {y, 2}] - 2*b),
D[c, t] + I*f[y]*c + I*d == 1/m*(D[c, {y, 2}] - 2*c)};
ic = {a == a0[y], b == b0[y], c == c0[y], d == d0[y]} /. t -> 0;
bc = {{a == 0, b == 0, c == 0, d == 0} /.
y -> -1, {a == 0, b == 0, c == 0, d == 0} /. y -> 1};];
ptoafunc = pdetoae[{a, b, c, d}[t, y], grid /@ {t, y}, difforder];
del = #[[2 ;; -2]] &;
ae = Map[del, Most /@ ptoafunc[eq], {2}];
aeic = Map[del, ptoafunc[ic]];
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {a, b, c, d}, grid[t], grid[y]];
{barray, marray} =
CoefficientArrays[Flatten[{ae, aeic, aebc}], Flatten[var]];
sollst = LinearSolve[marray, -barray];
solmatlst = ArrayReshape[sollst, var // Dimensions];
solfunclst = ListInterpolation[#, grid /@ {t, y}] & /@ solmatlst;


Update:

pdetoae can solve this system, however, the solution obtained is totally noisy, like the following figure. Also I found if using ic = {a == a0[y], b == b0[y], c == c0[y], d == d0[y]} /. t -> Tend, that is, using Tend as the initial time, the solution seems reasonable. I am quite sure that the problem stems from some mistakes in the removal of boundary difference equations for implementation of BCs/ICs instead of using higher difference order. According to this answer, I understood that: the b.c.s of the problem are Dirichlet type so there's no need to use pdetoae on them, as did in the code.

Some specific questions:

1. In ae = Map[del, Most /@ ptoafunc[eq], {2}], why we should apply del to each element on the 2nd level in list fptoafunc[feq] with its last element removed (by Most[list])?

2. How to store the solution, which can then be used to solve for another equation? For example,

NDSolve[{D[u[t, y], t] - 1/m*D[u[t, y], {y, 2}] == -2*Real[b[t, y]*Conjugate[D[a[t, y], y]] + I*Conjugate[c[t, y]]*a[t, y]], u[0, y] == 0, u[t, -1] == 0, u[t, 1] == 0}, u, {t, 0, Tend}, {y, -1, 1}]

in which the solution a[t,y], b[t, y] and c[t, y] should be used to compute the RHS of the diffusion equation.

Could anybody help me? Thank you for any suggestion.

• 1. Please read this comment carefully: mathematica.stackexchange.com/questions/174772/… Notice pdetoae is merely auxiliary function for the implementation of FDM, you must first understand the basic of FDM before using pdetoae. 2. "I didn't want to use Dumpsave since it normally is used to save the full info. " What do you mean by "full info"? – xzczd Jan 22 '19 at 11:10
• @xzczd 1. already read all the comments again and again but cannot figure out mistakes. I don't understand Most /@ ptoafunc[eq] although I knew the basic FDM. 2. By ''full info.'' I meant all data of an InterpolationFunction, including mesh, all the derivative values, etc. Here, only the time evolution of the original functions and the first space-derivatives are needed. If my understanding is wrong, please kindly suggest. Thanks in advance. – user55777 Jan 22 '19 at 11:49
• 1. Think about how to solve $\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2},u(0,x)=x(1-x),u(t,0)=u(t,1)=0$ with pdetoae first. 2. Your guess for InterpolatingFunction is wrong. Just check Interpolation[{1, 3, 6, 4, 5}] // InputForm. Also, check this post. – xzczd Jan 22 '19 at 12:07
• "Since it is an explicit method, having nothing to do with LinearSolve. " No, even if it's an explicit scheme, one can still write down every equation and solve the whole system all in once with LinearSolve. Yes, the answer for your question is replace Most with Rest, now tell me why. – xzczd Jan 23 '19 at 11:23
• As I've already mentioned, InterpolatingFunction isn't as expensive as you've imagined. Store it with DumpSave is the proper way to go. – xzczd Jan 23 '19 at 13:06

## 1 Answer

All issues are resolved by @xzczd. I will publish the solution to the problem so that the topic is completed

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,
MapThread[
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]]

Clear[f, m, Tend];
f[y_] := y; m = 300; Tend = 10;

(*Ramdom Iintial Conditions*)
Clear[a0, b0, c0, d0]

SeedRandom;
a0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
b0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
c0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];
SeedRandom;
d0 = Interpolation[
Table[{(-20 + i)/20, If[1 <= i < 40, RandomReal[{-1, 1}], 0]}, {i,
0, 40}], InterpolationOrder -> 4];

Clear[difforder, domain, points, grid];

difforder = 4;
domain[y] = {-1, 1}; domain[t] = {0, Tend};
points[y] = 201; points[t] = 101;
(grid@# = Array[# &, points@#, domain@#]) & /@ {y, t};

Clear[a, b, c, d];
With[{a = a[t, y], b = b[t, y], c = c[t, y], d = d[t, y]},
eq = {I*a + D[b, y] + I*c == 0,
D[a, t] + I*f[y]*a + b*D[f[y], y] + I*d ==
1/m*(D[a, {y, 2}] - 2*a),
D[b, t] + I*f[y]*b + D[d, y] == 1/m*(D[b, {y, 2}] - 2*b),
D[c, t] + I*f[y]*c + I*d == 1/m*(D[c, {y, 2}] - 2*c)};
ic = {a == a0[y], b == b0[y], c == c0[y], d == d0[y]} /. t -> 0;
bc = {{a == 0, b == 0, c == 0, d == 0} /.
y -> -1, {a == 0, b == 0, c == 0, d == 0} /. y -> 1};];
ptoafunc = pdetoae[{a, b, c, d}[t, y], grid /@ {t, y}, difforder];
del = #[[2 ;; -2]] &;
ae = Map[del, Rest /@ ptoafunc[eq], {2}];
aeic = Map[del, ptoafunc[ic]];
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {a, b, c, d}, grid[t], grid[y]];
{barray, marray} =
CoefficientArrays[Flatten[{ae, aeic, aebc}], Flatten[var]];
sollst = LinearSolve[marray, -barray];
solmatlst = ArrayReshape[sollst, var // Dimensions];

solfunclst = ListInterpolation[#, grid /@ {t, y}] & /@ solmatlst

sol = NDSolveValue[{D[u[t, y], t] - 1/m*D[u[t, y], {y, 2}] == -2*
Re[solfunclst[][t, y]*
Conjugate[Derivative[0, 1][solfunclst[]][t, y]] +
I*Conjugate[solfunclst[][t, y]]*solfunclst[][t, y]],
u[0, y] == 0, u[t, -1] == 0, u[t, 1] == 0},
u, {t, 0, Tend}, {y, -1, 1}];
p = {"a", "b", "c", "d"};

{Table[Plot[ReIm@solfunclst[[i]][10, y], {y, -1, 1},
PlotLabel -> p[[i]]], {i, 4}],
Plot3D[sol[t, y], {t, 0, Tend}, {y, -1, 1}, Mesh -> None,
ColorFunction -> Hue, AxesLabel -> {"t", "y", "u"}]}
` 