Solve a nonlinear PDE equation with a Neumann boundary condition

Question

I am trying to use Mathematica 10 to solve a PDE $$u_t=u_{xx}+u_{yy}+u(1-u),$$ in the unit disk $(x,y) \in D=\{(x,y):x^2+y^2<1\}$, with the Neumann boundary condtion $$\frac{\partial u}{\partial n} \bigg|_{\partial D}=0,$$ and the initial condition $$u(x,y,0)=x^2+y^2-1.$$ My Mathematica code is the following

Ω = ImplicitRegion[x^2 + y^2 <= 1, {x, y}];
RegionPlot[Ω, AspectRatio -> Automatic]

Subscript[Γ, N] = NeumannValue[0, x^2 + y^2 == 1];
op = 
  D[u[x, y, t], t] - Laplacian[u[x, y, t], {x, y}] - u[x, y, t] (2 - u[x, y, t]);
uif = 
  NDSolveValue[
    {op == Subscript[Γ, N],u[x, y, 0] == x^2 + y^2 - 1}, 
    u, {x, y} ∈ Ω, {t, 0, 10}]

But when I evaluate the code, Mathematica gives the error message

Nonlinear coefficients are not supported in this version of NDSolve.

What do I need to do? Any help is appreciated.

Answer 1

This problem can (mostly) be solved if we convert to polar coordinates and solve u as a function of r,t. The pde in the narrative does not match the code, but I will go with the code, changing the variables from x,y to r.

pde = Laplacian[u[r, t], {r, theta]}, "Polar"] + u[r, t] (2 - u[r, t]) == D[u[r, t], t]
(*D[ u[r, t], r]/r + D[ u[r, t], r, r] + (2 - u[r, t]) u[r, t] == 
 D[ u[r, t], t]*)

NDSolve doesn't like the 1/r in the pde at r = 0, so instead of solving from 0, I will have r go from some small value epsilon to 1. In doing that, I will have to add the boundary condition at r = epsilon. The derivative of u would be 0 at the center of the disk anyway.

epsilon = .0001;
bc1 = (D[ u[r, t], r] /. r -> 1) == 0;
bc2 = (D[ u[r, t], r] /. r -> epsilon) == 0;

The given initial condition does not match up with the boundary condition at r = 1 very well. NDSolve in trying to reconcile the two conditions is too unstable to get anywhere, so I modified the ic slightly, by making it a combination of 2 functions. Most of the disk retains the given ic, with a transition to a funtion that has a 0 derivative at r = 1.

f1 = r^2 - 1
f2 = -1000 (r - 1)^2

Find the transition point, where the two are equal:

r1 = r /. FindRoot[f1 - f2, {r, .9}]
(*0.998002*)

Form the ic transition from f1 to f2 with a combination of UnitStep functions.

ic = u[r, 0] == (UnitStep[r - epsilon] - 
      UnitStep[r - r1]) f1 + (UnitStep[r - r1] - UnitStep[r - 1]) f2

Now plug into NDSolve.

s = NDSolve[{pde, bc1, bc2, ic}, u, {r, epsilon, 1}, {t, 0, 1}, 
  PrecisionGoal -> Infinity, MaxStepFraction -> 1/1000, 
  MaxSteps -> {50000, Automatic}];

NDSolve still complains about inconsistent bc's and ic's and doesn't run all that long, but it gives a enough of a solution to see the behavior of u.

u[r_, t_] = u[r, t] /. s[[1]];

gifs = Table[Plot[u[r, t], {r, epsilon, 1}, PlotRange -> {-2, 0}, 
    PlotLabel -> t "t"], {t, 0, .5, .01}];
ListAnimate[gifs]

Answer 2

Because of the differential boundary condition the circular eigenfunctions, which result from a separation of variables approach, do not have linear Fourier series coefficients. You will likely need to look for the cross points in the eigenvalues of two trigonometric functions. To figure out which ones assume that there is a function of U is explicitly a function of X and Y, say U(x,y)=G(x)H(y). Substitute into the PDE, then put all of the x variables and time on one side of the equation, then all of the y variables on the other side of the equation. Set each side of the differential equation equal to a constant, lambda, which is the eigenvalue. Solve the ordinary differential equations for G(x) and H(y) in terms of the eigenvalue separately. The differential equations will have different functional solutions. For a hint, you can look at a similar ODEs here: https://en.wikipedia.org/wiki/Fisher%27s_equation.
Plug the Fourier series solution of Sum lambda(n)*G(x)*H(y) into the boundary condition and the initial condition to eliminate some eigenfunctions and find the relationship between the eigenvalues. You'll need to find the eigenvalues within mathematica manually, then compose an infinite series of it and the eigenfunctions. Then plot the Fourier series on the disk using a Contour plot. You can typically expect to run into exponential or trigonometric eigenfunctions.