I am attempting to solve for waves on a water surface starting with a two dimensional solution. The equations are that the water must satisfy Laplace's equation everywhere with a time dependent boundary condition on the top surface. The equations are

$$\begin{align*} &\nabla^2\phi=0\\ &\left(\frac{\partial\phi}{\partial y}+\frac1{g}\frac{\partial^2\phi}{\partial t^2}\right)_{y=y_2}=0 \end{align*}$$

where $x$ is horizontal and $y$ vertical; also, $g$ is the acceleration due to gravity. The second equation applies at the surface. The first term in this equation is the vertical gradient so I think I can use a Neumann condition to enter this equation.

So I start by making a grid.

x2 = 4; y2 = 1;
reg = ImplicitRegion[0 <= x <= x2 && 0 <= y <= y2, {x, y}];
mesh = ToElementMesh[reg, "BoundaryMeshGenerator" -> {"Continuation"},
    MaxCellMeasure -> .002, "MaxBoundaryCellMeasure" -> 0.01];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All]

Mathematica graphics

I now make some initial conditions. I have worked out an initial velocity for the top surface and solve to give the initial value of the potential function.

a = 0.40825787026798765;
b = 0.1689925306573793;
g = 9.81;

solIC = NDSolveValue[{
    Laplacian[ϕ0[x, y], {x, y}] ==
     NeumannValue[-b (E^-(x - 1)^2 - a), 0 <= x <= x2 && y == y2],
    DirichletCondition[ϕ0[x, y] == 0, x == 0 && y == 0]
    }, ϕ0, {x, y} ∈ mesh];

The velocity on the surface and throughout the water looks reasonable.

Plot[Evaluate[(D[solIC[x, y], y]) /. y -> y2], {x, 0, x2}, 
 AspectRatio -> 1/4]
f[x_, y_] := Evaluate[Grad[solIC[x, y], {x, y}, "Cartesian"]];
StreamPlot[f[x, y], {x, 0, x2}, {y, 0, y2}, AspectRatio -> Automatic]

Mathematica graphics

Now I attempt to set up the time dependent problem.

sol = NDSolveValue[{
   Laplacian[ϕ[x, y, t], {x, y}] ==
    NeumannValue[1/g Derivative[0, 0, 2][ϕ][x, y, t], 
     0 <= x <= x2 && y == y2],
   DirichletCondition[ϕ[x, y, t] == 0, x == 0 && y == 0],
   ϕ[x, y, 0] == solIC[x, y]
   }, ϕ,
  {x, y} ∈ mesh, {t, 0, 0.1}

This does not work and gives me the error

CoefficientArrays::poly: ϕ$1664+ϕ$1665-NeumannValue[0.101937 ϕ$1666,0<=x<=4&&y==1] 
is not a polynomial. >>

Now I am stuck. This does not help. Any suggestions?

Suggestions from comments

  1. If I put in a time derivative in the equation this gives the same error.

  2. If I take the time derivative out of the NeumannValue and put a time derivative into the equation then it solves but is meaningless. Does this suggest we can't have time derivatives in the boundary conditions?

Can you suggest how we might set up an alternative equation who's solution might go into the boundary condition?

sol = NDSolveValue[{
    Laplacian[ϕ[x, y, t], {x, y}] == 0.1 D[ϕ[x, y, t], t] +
      NeumannValue[ϕ[x, y, t], 0 <= x <= x2 && y == y2],
    DirichletCondition[ϕ[x, y, t] == 0, x == 0 && y == 0],
    ϕ[x, y, 0] == solIC[x, y]
    }, ϕ,
   {x, y} ∈ mesh, {t, 0, 0.1}
Plot3D[sol[x, y2, t], {x, 0, x2}, {t, 0, 0.1}]

Mathematica graphics

So the problem could be the derivative in the NeumannValue.

  • $\begingroup$ I am not quite sure, but it seems that the equation is not of the type supported by Mma, just for the reason that the time derivative only enters the Neumann condition, but not the equation itself. May be one could try to add the first derivative multiplied by some factor to the equation as a regularization. And then try to get this factor down? $\endgroup$ Aug 3, 2016 at 13:31
  • $\begingroup$ @AlexeiBoulbitch I added in a time derivative to the equation and got the same error. So no progress. Thanks for the suggestion. $\endgroup$
    – Hugh
    Aug 3, 2016 at 14:02
  • $\begingroup$ Let's ignore the NeumannValue for a second, should there not be a Derivative[0, 0, 1][\[Phi]][x, y, t] or D[\[Phi][x, y, t], {t, 1}] in the time dependent equation? Also is this a wave equation? (second order time deriv? - but maybe the title is miss leading me here.) $\endgroup$
    – user21
    Aug 3, 2016 at 14:46
  • 1
    $\begingroup$ @user21 The text book is: Landau and Lifshitz, Course of Theoretical Physics, Fluid Mechanics 1987 page 32. They do analytical solutions. It is a linearised equation. I also found it here web.mit.edu/2.20/www/lectures/lec14/lecture14.pdf in section 6.2.1. $\endgroup$
    – Hugh
    Aug 3, 2016 at 16:34
  • 1
    $\begingroup$ @Hugh, I have an idea :-) $\endgroup$
    – user21
    Dec 2, 2016 at 12:08

1 Answer 1


Update: When I wrote this answer it was not clear to me that this question is about Airy wave theory which I do not know how to solve.

Old answer:

I am not sure if this is what you are looking for but here is a time dependent solution. We can take it from here if it's not what you need.

sol = NDSolveValue[{1/g D[ϕ[x, y, t], {t, 1}] - 
      Laplacian[ϕ[x, y, t], {x, y}] == 0, 
    DirichletCondition[ϕ[x, y, t] == 0, 
     x == 0 && y == 0], ϕ[x, y, 0] == solIC[x, y]}, ϕ, {x, 
     y} ∈ mesh, {t, 0, 0.1}];
Plot3D[Evaluate[{sol[x, y, 0], sol[x, y, 0.1]}], {x, y} ∈ 
  mesh, PlotLegends -> Automatic]

enter image description here

You can have derivatives of time in the NeumannValue - but that needs to be a wave type of equation:

sol = NDSolveValue[{1/g D[ϕ[x, y, t], {t, 2}] - 
      Laplacian[ϕ[x, y, t], {x, y}] == 
     NeumannValue[-Derivative[0, 0, 1][ϕ][x, y, t], 
      0 <= x <= x2 && y == y2], 
    DirichletCondition[ϕ[x, y, t] == 0, 
     x == 0 && y == 0], ϕ[x, y, 0] == solIC[x, y], 
    Derivative[0, 0, 1][ϕ][x, y, 0] == 0}, ϕ, {x, 
     y} ∈ mesh, {t, 0, 0.1}];

For this example, however, you'd need to set up an initial velocity that makes sense. Setting it to zero is not correct.

Here is a simplified example of how to set time derivatives within the NeumannValue. This is in essence an absorbing boundary condition (Note, that in this 1D case there is a small bug that will be fixed in a upcoming version, the sign in the NeumannValue is not correct. This has no effect on the 2D case you are looking at):

eqn = D[u[t, x], {t, 2}] == 
   D[u[t, x], {x, 2}] + 
     If[$VersionNumber <= 10.4, 1, -1]*Derivative[1, 0][u][t, x], 
     x == 0 || x == 1];
u0[x_] := Evaluate[D[0.125 Erf[(x - 0.5)/0.125], x]];
ic = {u[0, x] == u0[x], Derivative[1, 0][u][0, x] == 0};
ufun = NDSolveValue[{eqn, ic}, u, {t, 0, 1}, {x, 0, 1}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement"}}];
list = Table[
   Plot[ufun[t, x], {x, 0, 1}, PlotRange -> {-0.1, 1.3}], {t, 0, 1, 

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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