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.


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;
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.

  • $\begingroup$ Do you mean you want to define a custom method for NDSolve to solve the PDE system, similar to those introduced in tutorial/NDSolveExplicitRungeKutta? If so, I should say it's very hard, if not impossible. ExplicitRungeKutta seems to be the only method that can be highly customized in NDSolve, at least according to the introduction in the document. That's the reason I create pdetoode, fix etc. attempting to break the limitation somewhat outside of NDSolve. BTW, your problem is similar to this one: mathematica.stackexchange.com/q/163923/1871 The bottleneck is eq[[1]]. $\endgroup$
    – xzczd
    Oct 13, 2018 at 8:39
  • $\begingroup$ Thanks@xzczd Reply: R1. nope, I didn't mean to define a totally custom method. I want to explicitly discretize the spatial derivative with pdetoode and fix the order of the derivatives, then solve the obtained ODEs. R2. yes, eq[[1]] is hard to handle due to the absence of time derivative. From above link I realized that D on eq[[1]] allows NDSolve to choose an ODE solver. Questions: Q1. In pdetoode, can u explicitly choose the spatial scheme, say, 2nd-order central difference? Did u mean it's impossible or very hard to control time-marching scheme in NDSolve even with pdetoode? $\endgroup$
    – user55777
    Oct 13, 2018 at 11:13
  • $\begingroup$ 1. When the 4th parameter of pdetoode is set to 2, 2nd order central difference formula is chosen for all the grid points except for those on the boundary, you can play with a coarse grid, say, Array[#&, 5, {0,1}] and observe the resulting system. 2. By "customize" I mean creating something not built-in. If you just want to control the ODE solver, it's not hard, but do notice the ODE solver is robust enough and don't need any manual adjustion in most cases, the bottleneck for PDE solving is almost always the spatial discretization. $\endgroup$
    – xzczd
    Oct 13, 2018 at 11:38
  • $\begingroup$ Also, notice a simple D[eq[[1]], t] won't invoke ODE solver, because the system still doesn't involve D[q, t] term afterwards. $\endgroup$
    – xzczd
    Oct 13, 2018 at 12:05
  • 2
    $\begingroup$ What is the actual question you have? $\endgroup$
    – user21
    Oct 13, 2018 at 12:18

1 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, 

enter image description here

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.

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Jan 22, 2019 at 10:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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