11

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.

user55777
  • 671
  • 3
  • 16
  • 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: https://mathematica.stackexchange.com/q/163923/1871 The bottleneck is eq[[1]]. – xzczd Oct 13 '18 at 08:39
  • 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? – user55777 Oct 13 '18 at 11:13
  • 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.
    • – xzczd Oct 13 '18 at 11:38
  • 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. – xzczd Oct 13 '18 at 12:05
  • 2
    What is the actual question you have? – user21 Oct 13 '18 at 12:18
  • @xzczd as for Q2 in my post, to integrate PDE from endtime to 0, don't we need additional modifications in NDSolve or/and pdetoode, except modifying ic and time limit in {t, endtime, 0}? – user55777 Oct 13 '18 at 14:14
  • Well, if you have difficulty in understanding the source code of pdetoode, then it's not too bad an idea to let it alone. You just need to know it's usage. As to the removing step, notice ptoofunc generates discretized equation on every grid points, so we need to delete some of them to make room for those i.c.s/b.c.s. Is your basic goal solving this PDE system? If so, I can post an answer. – xzczd Oct 13 '18 at 14:17
  • @user21 the actual question is: for a system of pdes with given 'ICs' at endtime, how to find the functions at t=0 by explicitly encoding a FD scheme in NDSolve. – user55777 Oct 13 '18 at 14:19