Solving systems of partial differential equations with functions of different number of variables

Question

I am trying to solve the following system of two partial differential equations

$\partial_t G(x,y,t) + \partial_x G(x,y,t)+\partial_y G(x,y,t) = -i\left[f(x,t) + f(y,t)\right] $

$\partial_t f(x,t) + \partial_x f(x,t) = f(x,t) - iG(x,0,t)$

with initial conditions $f(x,0)=0$ and $G(x,y,0) = G_0(x,y)$ where $G_0$ is a given function (e.g., a double Gaussian packet in the x-y plane)

I tried the following code:

L = 1; T = 1; x0 = -L/4; sigma = L/30; 
sol = NDSolve[{
   D[G[x, y, t], t] == -(D[G[x, y, t], x] + D[G[x, y, t], y]) - I*(f[x, t] + f[y, t]),
    D[f[x, t], t] == -D[f[x, t], x] + f[x, t] - I*(G[x, 0, t]),
    G[x, y, 0] == Exp[-((x - x0)/(Sqrt[2]*sigma))^2 - ((y - x0)/(Sqrt[2]*sigma))^2],
    f[x, 0] == 0
   },
  {G, f}, {x, -L/2, L/2}, {y, -L/2, L/2}, {t, 0, T},
  MaxSteps -> 500]`

But I get the errors:

Function::fpct: Too many parameters in {x,y,t} to be filled from Function[{x,y,t},0][x,t].

Function::fpct: Too many parameters in {x,y,t} to be filled from Function[{x,y,t},0][-0.5,0.].

Function::fpct: Too many parameters in {x,y,t} to be filled from Function[{x,y,t},0][-0.5,0.].

General::stop: Further output of Function::fpct will be suppressed during this calculation.

NDSolve::ndnum: Encountered non-numerical value for a derivative at y == 0.`.

NDSolve::ndnum: Encountered non-numerical value for a derivative at y == 0.`.

Any suggestion on how to proceed?

Answer 1

First of all, I'd like to point out the warning Function::fpct is probably because you've executed f = Function[{x, y, t}, 0]; and forgotten to Clear[f] then. This isn't the main issue, of course.

The main issue is, currently NDSolve can't handle equation system whose unknown functions are defined on different domains. (In your case, $G\in [-\frac{L}{2},\frac{L}{2}]\times[-\frac{L}{2},\frac{L}{2}]\times[0,T]$ and $f\in [-\frac{L}{2},\frac{L}{2}]\times[0,T]$. ) So let's discretize the system to an ODE system ourselves.

I'll use pdetoode for the generation of ODEs.

First, supplement b.c. to the system. Since you've mentioned b.c. isn't important here, I simply use zero Dirichlet b.c.:

L = 1; T = 1; x0 = -L/4; sigma = L/30;
domain = {-L/2, L/2};
{eq1, eq2} = {D[G[x, y, t], t] == -(D[G[x, y, t], x] + D[G[x, y, t], y]) - 
     I (f[x, y, t] + f2[x, y, t]), 
   D[f[x, t], t] == -D[f[x, t], x] + f[x, t] - I (G2[x, t])};

{ic1, ic2} = {G[x, y, 0] == 
    Exp[-((x - x0)/(Sqrt[2] sigma))^2 - ((y - x0)/(Sqrt[2] sigma))^2], f[x, 0] == 0};

{bc1, bc2} = {G[x, y, t] == 0 /. Outer[{# -> #2} &, {x, y}, domain],
   f[x, t] == 0 /. List /@ Thread[x -> domain]};

Notice I've modified the form of PDE system a bit (f[x, t] -> f[x, y, t], f[y, t] -> f2[x, y, t], G[x, 0, t] -> G2[x, t] ) because pdetoode can't handle functions defined in different domains all in once, either.

Next step is discretization. I've defined 2 functions ptoofunc1 and ptoofunc2 for the discretization of 2 domains.

points = 71; 
grid = Array[# &, points, domain];
difforder = 4;
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)

ptoofunc1 = pdetoode[{G, f, f2}[x, y, t], t, {grid, grid}, difforder];
ptoofunc2 = pdetoode[{G2, f}[x, t], t, grid, difforder];

del = #[[2 ;; -2]] &;
rule1 = {f[x_, y_][t_] :> f[x][t], f2[x_, y_][t_] :> f[y][t]};
rule2 = G2[x_][t_] :> G[x, 0][t];
ode1 = del /@ del@ptoofunc1@eq1 /. rule1;
ode2 = del@ptoofunc2@eq2 /. rule2;
odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;

Notice points should be an odd number, or G[x, 0, t] won't be properly handled.

del is a function for deleting equations at boundary to "make room" for b.c.s because ptoofunc1 and ptoofunc2 generate equations for every grid points. rule1 and rule2 is for transforming f2 and G2 back to f and g.

diff = With[{sf = 1}, D[#, t] + #] &;
odebc1 = Map[diff, MapAt[del /@ # &, ptoofunc1@bc1, {1}], {-2}];
odebc2 = Map[diff, ptoofunc2@bc2, {-2}];

diff is for transforming the disretized b.c. to (almost) equivalent ODEs. Well I admit the code above is a bit advanced, to have a better understanding for the whole process you may want to read this post.

The last step is to solve the system and rebuild the solutions to 2 interpolating functions:

sol = NDSolveValue[{ode1, ode2, odeic1, odeic2, odebc1, odebc2}, {Outer[G, grid, grid], 
    f /@ grid}, {t, 0, T}];

solG = rebuild[sol[[1]], {grid, grid}, 3];
solf = rebuild[sol[[2]], grid, 2];

Let's check the solution:

Manipulate[Plot3D[solG[x, y, t] // Evaluate , {x, ##}, {y, ##}, 
        PlotRange -> {-0.1, 1}], {t, 0, T}] & @@ domain

Manipulate[Plot[solf[x, t] // Abs // Evaluate, {x, ##}, PlotRange -> {0, 0.2}], {t, 0, 
    T}] & @@ domain

If you prefer periodic b.c.:

ptoofunc1 = pdetoode[{G, f, f2}[x, y, t], t, {grid, grid}, difforder, True];
ptoofunc2 = pdetoode[{G2, f}[x, t], t, grid, difforder, True];

ode1 = ptoofunc1@eq1 /. rule1;
ode2 = ptoofunc2@eq2 /. rule2;
odeic1 = ptoofunc1@ic1;
odeic2 = ptoofunc2@ic2;

sol = NDSolveValue[{ode1, ode2, odeic1, odeic2}, {Outer[G, grid, grid], f /@ grid}, {t, 
    0, T}];

solG = rebuild[sol[[1]], {grid, grid}, 3];
solf = rebuild[sol[[2]], grid, 2];

Answer 2

Well from the link below DSolve , in this stage of development can only handle 2 independent variables. http://reference.wolfram.com/language/tutorial/DSolveIntroductionToPDEs.html