2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – xzczd Nov 16 '17 at 4:26
  • $\begingroup$ 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. $\endgroup$ – Michele Cotrufo Nov 16 '17 at 16:26
2
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here enter image description here

$\endgroup$
0
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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