2nd order PDE: NDSolve error "Boundary condition is not specified on a single edge of the boundary of the computational domain"

Question

I am trying to solve a second order PDE. From the physics point of view all the boundary and Neumann conditions listed are correct, but still I get errors. Here is the code:

NDSolve[{I D[u[z, x, y], z]+D[u[z, x, y], x, x]+D[u[z, x, y], 
y,y]-(7*u[z,x,y])/(1+Abs[u[z,x,y]]^2) ==0 , u[0, x, y] == E^-(x^2 + 
y^2),(u^(1,0,0))[z, 0, 0] == 0}, u, {z, 0, 4}, {x,-5, 5}, {y, -5, 
5},Method->{PDEDiscretization->FiniteElement}]

If necessary, the following boundary conditions can be added also, but it doesn't really help.

u[z, -5, y] == 0, u[z, 5, y] == 0, u[z, x, -5] == 0, u[z, x, 5] == 0

Can anyone help to deal with this issue?

Welcome! I suggest the following:

As you receive help, try to give it too, by answering questions in your area of expertise.

Take the tour and check the faqs!

When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! — , Aug 24 '16 at 06:50
See How to copy code from Mathematica so it looks good on this site and Hints and help on using the site efficiently and How do I format my posts using Markdown or HTML? and read the FAQ — , Aug 24 '16 at 06:52

xzczd · Accepted Answer · 2022-06-23T06:37:34.147

Solution for user before v11 (v10?)

So this is another clue that NDSolve is silently improved (bug-fixed?) through the years. user21's solution works in v11 (tested on cloud) but does fail in v9.0.1 with warning ndcf and eerr generated (BTW in v8.0.4 the calculation simply never ends), and I can't find a way to make NDSolve solve it directly, so let's discretize the PDE to a set of ODEs with pdetoode (its definition can be found here).

The apparently erroneous condition (u^(1,0,0))[z, 0, 0] == 0 has been taken away.

eqn = I D[u[z, x, y], z] + D[u[z, x, y], x, x] + D[u[z, x, y], y, y] - (7 u[z, x, y])/(
    1 + Abs[u[z, x, y]]^2) == 0;
ic = u[0, x, y] == E^-(x^2 + y^2);
bc = {u[z, -5, y] == 0, u[z, 5, y] == 0, u[z, x, -5] == 0, u[z, x, 5] == 0};
zend = 4;
xdomain = ydomain = {-5, 5};
xpoints = ypoints = 25;
xgrid = ygrid = Array[# &, xpoints, ydomain];
difforder = 4;
(* Definition of pdetoode isn't included in this code piece,
   please find it in the link above. *)
ptoo = pdetoode[u[z, x, y], z, {xgrid, ygrid}, difforder];
del = Delete[#, {{1}, {-1}}] &;
ode = del /@ del@ptoo@eqn;
odeic =(del/@del@)ptoo@ic;
odebc = With[{sf = 1}, MapAt[del, ptoo@diffbc[z, sf]@bc, {{1}, {2}}]];
var = Outer[u, xgrid, ygrid];
sollst = NDSolveValue[{ode, odeic, odebc}, var, {z, 0, zend}];
sol = rebuild[sollst, {xgrid, ygrid}];
Animate[Row[
  Plot3D[#@sol[z, x, y], {x, -5, 5}, {y, -5, 5}, PlotRange -> 1/2, 
     PerformanceGoal -> "Quality", ImageSize -> Medium] & /@ {Re, Im}], {z, 0, 4}]

Notice the "strange" definition for odebc is necessay, or NDSolveValue will spit out mconly and icfail and fail. To know more about this boundary condtion, check this obscure tutorial. (Particularly the part about Boundary Conditions. )

Thank you a lot for your reply. Which version of Mathematica do you use? I tried on v9.0 but it gives
MapAt::partw: Part{2} of pdetoode[...] does not exist NDSolveValue::dsvar: "!({z, 0, 4}["ValuesOnGrid"]) cannot be used as a variable" NDSolveValue::argm: NDSolveValue called with 0 arguments; 3 or more arguments are expected etc. — Cyd, Dec 04 '16 at 08:37
@cyd The function definition of pdetoode isn't included in this post, please find it in the link. — xzczd, Dec 04 '16 at 08:42

score 1 · Answer 2 · answered Aug 24 '16 at 13:16

This works:

NDSolveValue[{I D[u[z, x, y], z] + D[u[z, x, y], x, x] + 
    D[u[z, x, y], y, y] - (7*u[z, x, y])/(1 + Abs[u[z, x, y]]^2) == 0,
   u[0, x, y] == E^-(x^2 + y^2),
  u[z, -5, y] == 0, u[z, 5, y] == 0, u[z, x, -5] == 0, 
  u[z, x, 5] == 0
  }, u, {z, 0, 4}, {x, -5, 5}, {y, -5, 5}, 
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> "TensorProductGrid"}]

Note that I removed the derivative on z in the initial condition. Also I treated this as a time dependent problem. Since this a nonlinear equation the FEM in V11 can not handle that and you need to use the "TensorProductGrid". I have included the Dirichlet boundary conditions, you could also try something like this:

NDSolveValue[{I D[u[z, x, y], z] + D[u[z, x, y], x, x] + 
    D[u[z, x, y], y, y] - (7*u[z, x, y])/(1 + Abs[u[z, x, y]]^2) == 0,
   u[0, x, y] == E^-(x^2 + y^2),
  Derivative[0, 1, 0][u][z, -5, y] == 0,
  Derivative[0, 1, 0][u][z, 5, y] == 0,
  Derivative[0, 0, 1] u[z, x, -5] == 0,
  Derivative[0, 0, 1] u[z, x, 5] == 0
  }, u, {z, 0, 4}, {x, -5, 5}, {y, -5, 5}, 
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid"}}]

But it correctly will give a message that the initial condition does not match the boundary condition, but perhaps that's OK in your case as NDSolve will try to find initial conditions that match.

Tried the first option, but still there are errors. It says "NDSolveValue::ndcf: Repeated convergence test failure at z == 0.16946625650097352; unable to continue." and "Warning: scaled local spatial error estimate of 511.89139378820755 at z = 0.16946625650097352 in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 49 points may be too large to achieve the desired accuracy or precision. Etc." Could this be due to different versions of Mathematica? — Cyd, Aug 27 '16 at 14:53
The conditions on the second solution are not correct for my problem and still there are same errors displayed for another value of z. — Cyd, Aug 27 '16 at 14:58

2nd order PDE: NDSolve error "Boundary condition is not specified on a single edge of the boundary of the computational domain"

2 Answers2

Solution for user before v11 (v10?)