Nonlinear parabolic equations: finite difference method

[mercurial]

To the moderator: please move this to the section on differential equations if you think it would be better there.

I'm looking at a reaction-diffusion equation of the form

\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + f(u),

where, e.g., f(u) = u(1-u) for the Fisher-Kolmogorov eqation, and I'm assuming the spatial domain is the unit interval with Neumann boundary conditions:

\frac{\partial u}{\partial x}(0, t) = 0, \; \frac{\partial u}{\partial x}(1, t) = 0, \; t\geq 0.

To keep things totally simple, I'm using the backward Euler method. When the original equation is just the heat equation, i.e., f = 0, this method just involves solving a simple matrix equation,

BU^{n+1} = U^n,

for evolving the solution from time step n to n+1. The matrix of course depends on the boundary conditions. The nonlinear version of Euler's method involves a nonlinear multi-variable equation:

BU^{n+1} + F(U^{n+1}) = U^n,

which can be solved for each time step using Newton's method.

Where I'm getting confused is: exactly how do I incorporate the boundary conditions when deriving the system of nonlinear equations? My reasoning tells me that they should be the same as for the linear equation, i.e., the matrices B are the same for both the linear and nonlinear systems.

My solver works perfectly for the heat equation with Neumann conditions, but gives only garbage for the Fisher-Kolmogorov equation after the first time step (if it converges at all). I'm convinced that the problem is not a bug in programming Newton's method. Anyway I'm using a library for that, and I've checked the functions I've supplied a thousand times already. Anyone here have experience with nonlinear parabolic equations? They're not the kind of thing that's typically covered in basic numerical analysis texts.

 

[Emreth]

Can you show your equations you used for the boundary nodes?Also for inner nodes?i can check them for you.
How i do it is, use central diff to discretize bc's. You solve for the imaginary nodes (for 0, its one right outside say 0-, for 1 its 1+). Then you write the discrete form of the PDE then plug in the values for 1+ in terms of 1- and 0- in terms of 0+ (they are equal for your bcs). I used central diff because its 2nd order accurate in x, like the discrete form for the second derivative.I don't know if f(u) creates some weirdness but i don't think so, its a pretty well behaved function.
Also i use the link below, they have analytical solutions for many types of pde's. You can use that to verify your calculations.

http://eqworld.ipmnet.ru/index.htm

 

[mercurial]

I'm using FD equations at the boundary nodes exactly like you describe: "central difference" because it's 2nd order accurate. The discretization matrix looks like:

A =<br /> \begin{pmatrix}<br /> -2 &amp; 2 &amp; 0 &amp; \ldots \\<br /> 1 &amp; -2 &amp; 1 &amp; 0 &amp; \ldots \\<br /> 0 &amp; 1 &amp; -2 &amp; 1 &amp;\ldots<br /> \end{pmatrix},

and I assume you know the rest. Writing U^n = (U_0^n,\ldots,U_m^n)^T for the approximation at the n-time step, the backward Euler method for the heat equation is:

\frac{U^{n+1} - U^n}{k} = h^{-2}AU^{n+1},

which becomes

(I - (k/h^2)A)U^{n+1} = U^n.

The matrix I - (k/h^2)A is what I'm calling B. According to my reasoning, the nonlinear equation involves a function F:\mathbf{R}^{m+1}\to \mathbf{R}^{m+1} yielding a nonlinear backward Euler method:

\frac{U^{n+1} - U^n}{k} = h^{-2}AU^{n+1} + F(U^{n+1}).

In the case of the Fisher-Kolmogorov eqn., F_i(U) = U_i(1-U_i) in each coordinate. Again, I haven't found this in any books, but it's the most obvious thing to do. The system of equations I have to solve in each time step is:

BU^{n+1} -kF(U^{n+1}) - U^n = 0;

sorry for changing my notation from my first post a little. I have to supply the left hand side together with its Jacobian to the nonlinear solver. Here the "current" time step U^n is treated as a constant. The Jacobian of the LHS is B - kF&#039;, where the Jacobian F&#039; is the obvious diagonal matrix.

So far so good?

For Newton's method I'm using the multi-root solver in the GNU Scientific Library http://www.gnu.org/software/gsl/. For my initial condition I'm using a Gaussian,

u_0(r) = \exp(-(r-0.5)^2/\sigma^2),

with \sigma chosen so that the function is basically flat at the endpoints x_0 = 0 and x_m = 1. The solution I'm getting for the first time step clearly violates the boundary conditions, and taking this U^1 as a starting point, the solver doesn't even converge for the second time step.

The logical explanation is that I've done something wrong with the nonlinear solver, and I'll be embarrassed for wasting everyone's time if that turns out to be the case. But I've reexamined that code so often by now, I'm starting to entertain the idea that I've missed something subtle with the boundary conditions for the FD equations.

Thanks for your help in this matter!

Jerry

 

[Emreth]

It seems fine to me. one issue with Newton's method is you should be close enough to the root, otherwise it might not converge. maybe you might change the IC, if its ok, to something that satisfies the BC's also. then the guess values(the prev values) might be closer to the solution.

 

[Emreth]

Ah one more thing.first, you can evaluate the nonlinear terms in the current time step, solve the system with regular gauss elimination. Then use those values as the initial guess for the solution of the original equations where the nonlinear terms are evald in the forward time.

 

[mercurial]

Hi Emreth,

I said in my last post that the initial condition u_0 is "basically flat at the endpoints." Thinking that that may be the problem, I modified it to satisfy the Neumann BC's exactly: to no avail. The same problem persists.

My starting guess for finding U^{n+1} with Newton's method is (the logical choice) U^n itself. The time step k I'm using is so small that the evolution should be tame. I know how solutions of the Fisher-Kolmogorov equation behave, and behavior of my numerical solution on the interior of the domain agrees with that. The problem is just at the boundary, most notably at x_0 = 0.

Thanks too for the tip on solving the linear part of the system forward in time to get a better initial guess for the nonlinear equations. I'll take a look at that later and get back to you.

Jerry

 

[mercurial]

problem evidently resolved

It looks like I've found the problem, indeed a bug: a typo among variable names in initializing the system. Ludicrous how something like that can "emulate" a numerical error, especially ludicrous when considering how much time I spent looking for it! Thanks anyway to everyone who at least skimmed through the thread or thought about the problem.

Jerry