# Solving linear coupled PDEs by FDM

I am trying to solve some linear, coupled PDEs for perturbative analysis (first order in time, 3rd order in space), for which I then plan to take the global spatial maxima of their magnitudes and plot them across time to show the temporal evolutions of the individual perturbations.

To solve the equations, I have attempted to use the template provided by user xzczd.

p = .011;
ky = 10;
c = 27;
Ω = 2800/p;
L = p/0.3;
v = p^2* Ω/L;
A0 = 1;
xf=3;
Θ[x_] := -(3 c (1 + p^2 ky^2)/(2 p^2*v*ky)) (Sech[
x c/(2 p*Sqrt[c^2 - ky^2 A0^2])]^2);

With[{A = A[t, x], Θ1 = Θ1[t, x],
vy = vy[t, x]},

pde1 = D[
A + p^2 (ky^2 A +
2 A0*Θ[x]*D[Θ1, x] - D[A, x, x]),
t] - c*D[
A + p^2 (ky^2 A +
2 A0*Θ[x]*D[Θ1, x] -
D[A, x, x]), x] +
p^2 (Θ[x])^2 (D[A, t] -
c*D[A, x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) \
(2*(Θ[x])*D[A, x] + A0*D[Θ1, x, x]);

pde2 = A0 (D[Θ1, t] - c*D[Θ1, x]) -
p^2 (D[Θ[x], x] (D[A, t] - c*D[A, x]) +
D[A0 D[Θ1, x, x] + 2 Θ[x] D[A, x],
t] - c*D[
A0 D[Θ1, x, x] + 2 Θ[x] D[A, x],
x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) (2 A0 Θ[
x] D[Θ1, x] - D[A, x, x]) -
p^2 ky A0 D[vy, x, x];

pde3 = (1/(p^4 Ω^2)) (D[vy, t] - c*D[vy, x]) ==
ky A0^2 D[Θ1, x, x] +
D[2 ky A0 A Θ[x], x];]

pde = {pde1, pde2, pde3};

ic = {A[0, x] == 10^(-5),
vy[0, x] == 0, Θ1[0, x] == 0};

bc = {A[t, xf] == 0, (D[A[t, x], x] /. x -> xf) ==
0, (D[A[t, x], x] /. x -> 0) ==
0, (D[Θ1[t, x], x] /. x -> xf) ==
0, (D[Θ1[t, x], x] /. x -> 0) ==
0, Θ1[t, xf] == 0,
vy[t, xf] == 0, (D[vy[t, x], x] /. x -> 0) == 0};

begintime = 0; endtime = 80;
points@x = 200; points@t = 25;
m = 200;
difforder = 8;
domain@x = {0, 3}; domain@t = {begintime, endtime};
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, t};

ptoafunc =
pdetoae[{A, Θ1, vy}[t, x], grid /@ {t, x},
difforder];
del = #[[2 ;; -2]] &;
ae = Map[del, Most /@ ptoafunc@pde, {2}];
aeic = del /@ ptoafunc@ic;
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {A, Θ1, vy}, grid@t, grid@x];
{barray, marray} =
CoefficientArrays[Flatten@{ae, aeic, aebc},
Flatten@var]; // AbsoluteTiming
Block[{p = .011,
ky = 10,
c = 27,
Ω = 2800/p,
L = p/0.3,
v = p^2* Ω/L,
A0 = 1},
sollst = LinearSolve[N@marray, -barray]]; // AbsoluteTiming
solmatlst = ArrayReshape[sollst, var // Dimensions];
solfunclst = ListInterpolation[#, grid /@ {t, x}] & /@ solmatlst;
plot[f_, style_, t_] :=
Plot[solfunclst[t, x] // Through // f //
Evaluate, {x, domain@x} // Flatten // Evaluate,
PlotStyle -> style]


but I run into an error once I run the code:

{17.5364, Null}


Could someone help resolve the error in this code? Or if I should be approaching this problem through different means? I apologize for the general inexperience in this matter. Thank you very much in advance.

• Where's the definition for xf? – xzczd May 30 '19 at 11:01
• @xzczd My bad, I forgot to include it here. xf should equal 3. I still seem to encounter some errors, however.. – user213068 May 30 '19 at 11:05
• You've only given 5 b.c.s in the code (2 for A, 1 for vy, 2 for \[CapitalTheta]1), are you sure it's enough? Usually the number of b.c. should be equal to the highest differential order in the corresponding direction for every unknown function. – xzczd May 30 '19 at 11:13
• Notice usually the number of i.c./b.c. should be equal to the highest differential order in corresponding direction for every unknown function i.e. you need 3 i.c.s and 8 b.c.s in total. – xzczd May 30 '19 at 11:27
• As to the removing of redundant equations, since you have 3×25×200 unknown variables to solve, while 3×25×200 + 3×200 + 8×25 equations at hand, you need to remove 3×200 + 8×25 from the system. But are you sure now the system is correct? I've tried various settings for points, difforder, etc., but the solution always becomes unstable very fast. – xzczd May 30 '19 at 13:09

OK, let me extend my comments to an answer. Your code doesn't give proper result because you haven't removed the redundant equations properly.

First of all, notice pdetoae will discretize equations in the following way:

1. If the equation is defined on the whole domain of definition, difference equations will be generated on every grid points. In your case, you'll obtain 3 × points[t] × points[x] difference equations after discretizing the 3 PDEs.

2. If the equation is only defined on the boundary of the domain of definition, difference equations will only be generated on grid points on the boundary. In your case, you'll obtain 3 × points[x] equations after discretizing i.c.s, and 8 × points[t] equations after discretizing b.c.s.

Now here comes the problem, we only have 3 × points[t] × points[x] unknown variables in the discretized system, but 3 × points[t] × points[x] + 3 × points[x] + 8 × points[t] equations, so the system is over-determined and we need to remove some of the equations as redundant ones. Which ones should be removed? Those closest to the i.c.s and b.c.s. (But why? Honestly speaking, I don't know, but this strategy seems to always work well. )

The following is the fixed code. I've reduced endtime to 1/2 because the solution damps fast.

p = 0.011; ky = 10; c = 27; Ω = 2800/p; L = p/0.3; v = (p^2 Ω)/L; A0 = 1; xf = 3;
Θ[x_] := -(((3 c (1 + p^2 ky^2)) Sech[(x c)/(2 p Sqrt[c^2 - ky^2 A0^2])]^2)/(
2 p^2 v ky));
Clear[pde, ae, aeic]
With[{A = A[t, x], Θ1 = Θ1[t, x], vy = vy[t, x]},
pde@1 = D[A +
p^2 (ky^2 A + 2 A0*Θ[x]*D[Θ1, x] - D[A, x, x]), t] -
c*D[A + p^2 (ky^2 A + 2 A0*Θ[x]*D[Θ1, x] - D[A, x, x]),
x] + p^2 (Θ[x])^2 (D[A, t] -
c*D[A, x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) (2*(Θ[x])*D[A, x] +
A0*D[Θ1, x, x]);

pde@2 = A0 (D[Θ1, t] - c*D[Θ1, x]) -
p^2 (D[Θ[x], x] (D[A, t] - c*D[A, x]) +
D[A0 D[Θ1, x, x] + 2 Θ[x] D[A, x], t] -
c*D[A0 D[Θ1, x, x] + 2 Θ[x] D[A, x],
x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) (2 A0 Θ[
x] D[Θ1, x] - D[A, x, x]) - p^2 ky A0 D[vy, x, x];

pde@3 = (1/(p^4 Ω^2)) (D[vy, t] - c*D[vy, x]) ==
ky A0^2 D[Θ1, x, x] + D[2 ky A0 A Θ[x], x];
ic = {A == 1/10^5, Θ1 == 0, vy == 0} /. t -> 0;
bc = {{A == 0, D[A, x] == 0, D[Θ1, x] == 0, Θ1 == 0,
vy == 0} /. x -> xf, {D[A, x] == 0, D[Θ1, x] == 0,
D[vy, x] == 0}} /. x -> 0;]

begintime = 0; endtime = 1/2(*1/10*);
points@x = 200; points@t = 25;
difforder = 4;
domain@x = {0, xf}; domain@t = {begintime, endtime};
(grid@# = Array[# &, points@#, domain@#]) & /@ {x, t};

(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{A, Θ1, vy}[t, x], grid /@ {t, x}, difforder];
deletetwo = #[[2 ;; -2]] &;
deletethree = #[[2 ;; -3]] &;
{ae@1, ae@2} = deletethree /@ Rest@ptoafunc@pde@# & /@ {1, 2};
{ae@3} = deletetwo /@ Rest@ptoafunc@pde@# & /@ {3};
{aeic@1, aeic@2} = deletethree@ptoafunc@ic[[#]] & /@ {1, 2};
{aeic@3} = deletetwo@ptoafunc@ic[[#]] & /@ {3};
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {A, Θ1, vy}, grid@t, grid@x];
{barray, marray} =
CoefficientArrays[Flatten@{ae /@ {1, 2, 3}, aeic /@ {1, 2, 3}, aebc},
Flatten@var]; // AbsoluteTiming

sollst = LinearSolve[N@marray, -barray]; // AbsoluteTiming

solmatlst = ArrayReshape[sollst, var // Dimensions];
solfunclst = ListInterpolation[#, grid /@ {t, x}] & /@ solmatlst;

GraphicsRow[(Plot3D[#1[t, x], {t, begintime, endtime}, {x, 0, xf}, PlotRange -> All,
PlotPoints -> 50] &) /@ solfunclst] You may need to adjust difforder, points, etc. further to obtain an accurate enough result.

• Oh! That makes a lot of sense! Thank you so much, I really can't appreciate your help and patience enough! – user213068 May 30 '19 at 14:18
• @user213068 Glad it helps. Actually you don't need to accept that fast, it's OK to wait for a while to see if someone comes with a better answer. – xzczd May 30 '19 at 14:26

We can use the standard solver with special options. The pictures are not as beautiful as using pdetoae, but perhaps other solutions have been found here.

p = .011;
ky = 10;
c = 27;
\[CapitalOmega] = 2800/p;
L = p/0.3;
v = p^2*\[CapitalOmega]/L;
A0 = 1;
xf = 3;
\[CapitalTheta][
x_] := -(3 c (1 + p^2 ky^2)/(2 p^2*v*ky)) (Sech[
x c/(2 p*Sqrt[c^2 - ky^2 A0^2])]^2);

With[{A = A[t, x], \[CapitalTheta]1 = \[CapitalTheta]1[t, x],
vy = vy[t, x]},
pde1 = D[A +
p^2 (ky^2 A + 2 A0*\[CapitalTheta][x]*D[\[CapitalTheta]1, x] -
D[A, x, x]), t] -
c*D[A + p^2 (ky^2 A +
2 A0*\[CapitalTheta][x]*D[\[CapitalTheta]1, x] -
D[A, x, x]), x] +
p^2 (\[CapitalTheta][x])^2 (D[A, t] -
c*D[A, x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) \
(2*(\[CapitalTheta][x])*D[A, x] + A0*D[\[CapitalTheta]1, x, x]);
pde2 = A0 (D[\[CapitalTheta]1, t] - c*D[\[CapitalTheta]1, x]) -
p^2 (D[\[CapitalTheta][x], x] (D[A, t] - c*D[A, x]) +
D[A0 D[\[CapitalTheta]1, x, x] + 2 \[CapitalTheta][x] D[A, x],
t] -
c*D[A0 D[\[CapitalTheta]1, x, x] +
2 \[CapitalTheta][x] D[A, x],
x]) == ((v p^2 ky)/(1 + (p^2 ky^2))) (2 A0 \[CapitalTheta][
x] D[\[CapitalTheta]1, x] - D[A, x, x]) -
p^2 ky A0 D[vy, x, x];
pde3 = (1/(p^4 \[CapitalOmega]^2)) (D[vy, t] - c*D[vy, x]) ==
ky A0^2 D[\[CapitalTheta]1, x, x] +
D[2 ky A0 A \[CapitalTheta][x], x];]

pde = {pde1, pde2, pde3};

ic = {A[0, x] == 10^(-5), vy[0, x] == 0, \[CapitalTheta]1[0, x] == 0};

bc = {A[t, xf] == 0, (D[A[t, x], x] /. x -> xf) ==
0, (D[A[t, x], x] /. x -> 0) ==
0, (D[\[CapitalTheta]1[t, x], x] /. x -> xf) ==
0, (D[\[CapitalTheta]1[t, x], x] /. x -> 0) ==
0, \[CapitalTheta]1[t, xf] == 0,
vy[t, xf] == 0, (D[vy[t, x], x] /. x -> 0) == 0};

{as, vs, tets} =
NDSolveValue[{pde, ic, bc}, {A, vy, \[CapitalTheta]1}, {t, 0,
1}, {x, 0, xf},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 30, "MaxPoints" -> 30}}}];

{Plot3D[as[t, x], {t, 0, 1}, {x, 0, xf}, AxesLabel -> {"t", "x", ""},
PlotLabel -> "A", PlotRange -> All, Mesh -> None,
ColorFunction -> Hue],
Plot3D[vs[t, x], {t, 0, 1}, {x, 0, xf}, AxesLabel -> {"t", "x", ""},
PlotLabel -> "vy", PlotRange -> All, Mesh -> None,
ColorFunction -> Hue],
Plot3D[tets[t, x], {t, 0, 1}, {x, 0, xf},
AxesLabel -> {"t", "x", ""}, PlotLabel -> "\[CapitalTheta]1",
PlotRange -> All, Mesh -> None, ColorFunction -> Hue]} • Interesting. Just a side note: this solution won't work in v9.0.1, but works at least in v11.2 and v12. BTW, my solution isn't necessarily the better one, actually it varies quite a bit if points[t] becomes larger… – xzczd May 31 '19 at 17:39
• @xzczd Therefore, I have placed this code so that it can be seen that with the automatic selection of the integration step in t there may be other modes. – Alex Trounev May 31 '19 at 17:45
• @xzczd , Alex I actually had a question about the dependence of t with these methods of FDM. I was testing c=20.01, ky=20 and found that, while for large endtime (w. xzczd's algorithm), the solutions decayed very fast (endtime ~3), changing the endtime to zoom into this solution gives weird flaps for all the plots towards endtime (~0.2). I ran Alex's code for these parameters, but it returns a "scaled local spatial error estimate.. in the direction of independent variable x is much greater than the prescribed error tolerance." Is there a reason for these behaviours? Thanks so much – user213068 May 31 '19 at 19:54
• @AlexTrounev And thank you so much for your response! It's very helpful in learning the ways these methods may be employed. – user213068 May 31 '19 at 20:03
• @user213068 The eerr warning returns by Alex's solution indicates NDSolve thinks the obtained solution isn't accurate enough. As to the behavior of my solution, sadly stability analysis for certain finite difference scheme is beyond my (and most average user's, I'm afraid) reach and I don't have any theoretical explanation, judging whether the solution is reliable enough based on the physical background and by testing various parameters is the most practical way to go, I think. – xzczd Jun 1 '19 at 4:47