12
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ Greetings! Make the most of Mma.SE and take the tour now. Help us to help you, write an excellent question. Edit if improvable, show due diligence, give brief context, include minimal working examples of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$ Commented Feb 1, 2016 at 16:52
  • $\begingroup$ Unrelated to your problem, but you should avoid using Subscript while defining symbols (variables). Subscript[x, 1] is not a symbol, but a compound expression, if you do $x_1=2$ you are actually doing Set[Subscript[x, 1], 2] which is to assign a Downvalue to Subscript and not an Ownvalue to an indexed x as you may intend. Read how to properly define indexed variables here. $\endgroup$ Commented Feb 1, 2016 at 17:14
  • $\begingroup$ I realize that you give a simplified example, but if the problems you are interested in satisfy rotational symmetry you could try polar coordinates. $\endgroup$ Commented Aug 16, 2017 at 10:02
  • $\begingroup$ On a related note, in the subsection "Classical Partial Differential Equations" of this documentation page, it says in particular that the coefficient a (in your case a=1-u) is a scalar. I am guessing a can depend on {x,y}, but not on u, that's just not implemented in NDSolve. $\endgroup$ Commented Aug 20, 2017 at 15:04
  • $\begingroup$ Your PDE has the solution u=0. I guess that you are looking, however, for a nontrivial one, are you? Generally, MethodOfLines must work for this PDE. In the case of the equation like this one, it is not, however, straightforward. The possibility to apply this method depends on the solution you expect to get. $\endgroup$ Commented Oct 14, 2017 at 20:21

2 Answers 2

Reset to default
4
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
1
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for answering my question. I am seeking numerical solution for parAbolic PDEs with Neumann boundary conditions. My question is just a special example. For general parabolic PDEs, the method you mentioned is not working. $\endgroup$ Commented Feb 2, 2016 at 2:24

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.