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

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?

• The warning is probably because you forgot to Clear[f] first. Then you need to add proper artificial b.c. for your problem because otherwise we won't know what b.c. is added to NDSolve. (See this post for more information. ) The most troublesome part f[x, t] + f[y, t] still remains though, but I think it's not too hard to overcome. If you supplement proper b.c., I can have a try. – xzczd Nov 16 '17 at 4:26
• Hi, thanks for your answer. Using Clear[f] before the code did not remove the warning. Regarding the boundary conditions: I expect the initial "pulse" (in the xy plane) never to approach the boundaries, so probably the exact boundary conditions are not important. I think periodic boundary conditions should work, or whatever else you think it may be easy to implement. – Michele Cotrufo Nov 16 '17 at 16:26

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,

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];
`

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

• Thanks for the info. Do you know if the same limit applies to NDSolve? I could not find anything in the the documentation. – Michele Cotrufo Nov 16 '17 at 16:29