Solve PDEs with finite difference scheme by modifying NDSolve-based solver

Question

Motivation

As discussed here, NDSolve uses different difference orders for various spatial derivatives and the implicit design could cause trouble in certain cases. Thus, it may be useful to specify the difference order for spatial derivatives as well as customize difference scheme for time advance in some applications. And that's why I'm trying to use a finite difference method (FDM) encoded in NDSolve to construct a lower-level PDE solver instead of using the high-level NDSolve black box directly.

Problem

Consider a system of PDEs which will be solved with specified FDM encoded in myFDM[points, xdifforder, torder]:

secondDop[x_, α_] = (D[#, {x, 2}] - (α^2 + 1)*#) &; (*define a differential operator*)
With[{a = a[t, x], b = b[t, x], c = c[t, x], q = q[t, x], U = x},
 eq = {I*α*a - D[b, x] + I*c == 0,
   -D[a, t] == 1/m*secondDop[x, α][a] - I*α*U*a - I α q,
   -D[b, t] == -a*D[U, x] + 1/m*secondDop[x, α][b] - 
     I*α*U*b + D[q, x],
   -D[c, t] == 1/m*secondDop[x, α][c] - I*α*U*c - I*q};
 (* values in ICs could be adjusted to allow the code to be run *)
 ic = {a == 1-x^2, b == Sin[x π], c == 0, q == 1} /. t -> endtime;
 bc = {{a == 0, b == 0, c == 0, q == 1} /. 
    x -> -1, {a == 0, b == 0, c == 0, q == 1} /. x -> 1};]

endtime = 5; points = 100;
m=200;
xdifforder = 2; torder = 2; (* Difference orders in x and t *)
grid = Array[# &, points, {-1, 1}]; (*spatial points separated by a distance of 2/points*)
bsol=NDSolveValue[{eq, ic, bc}, {a, b, c, q}, {t, endtime, 0}, {x, -1, 1},
  Method -> [myFDM[points, xdifforder, torder]]]

Specifications and Questions

I believe the desired function myFDM should be a modification of the pdetoode solver developed here by xzczd with the following significant differences. (Although I have learned Wolfram language for a long time some parts in pdetoode are still confusing for me, e.g. the line to remove redundant equations from every PDE/IC in moveredundant...)

1. My problem includes 4 PDEs instead of one PDE;

2. These PDEs will be integrated from endtime to 0 due to mathematical reasons, i.e. backward, as implied by the minus sign before each time derivative term on the LHS.

Thank you very much.

Answer 1

As mentioned above, this problem is somewhat similar to, but even more troublesome than this one, because there seems to be no obvious way to discretize and transform the PDE system to an ODE system. To be more specific, it's hard to introduce derivative of $q$ with respect to $t$ without causing other trouble. So, let's give up NDSolve and turn to pure finite difference method (FDM). It's not too hard to achieve with the help of pdetoae.

Definitions for eq, ic, bc are the same as those in the question. α isn't specified in the question so I simply choose 1.

Clear[grid, points]
begintime = 0; endtime = 5; 
points@x = 200; points@t = 25;
m = 200;
difforder = 8;
domain@x = {-1, 1}; 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, b, c, q}[t, x], grid /@ {t, x}, difforder];
del = #[[2 ;; -2]] &;
ae = Map[del, Most /@ ptoafunc@eq, {2}];
aeic = del /@ ptoafunc@ic;
aebc = ptoafunc@bc;
var = Outer[#[#2, #3] &, {a, b, c, q}, grid@t, grid@x];
{barray, marray} = 
   CoefficientArrays[Flatten@{ae, aeic, aebc}, Flatten@var]; // AbsoluteTiming
Block[{α = 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,
   PlotRange -> {-2, 2}, PlotStyle -> style]
Manipulate[GraphicsRow@{plot[Re, Automatic, t], plot[Im, Dashed, t]}, {t, endtime, 
  begintime}]

Notice the setting difforder = 8 is important here, or the solution of q will be noisy.

Also, the behavior of q near the boundary seems to suggest the b.c.s for $q$ might be redudant.