Solving linear coupled PDEs by FDM
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:
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.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, and8 × 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,
please find it in the link above. *)
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.
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]}
