Linearization in Finite Difference Method: Why?

I have a very basic question, and I hope some of you might be able to help me: In Fluid Dynamics, a common equation turning up time and time again is the so-called continuity equation:

$$\frac{\partial}{\partial x}(T\frac{\partial h}{\partial x}) + R = S\frac{\partial h}{\partial t}$$

Where $T$ is the transmissivity of the soil matrix, defined by the product $Kh$, where $h$ is the hydraulic head and $K$ is the hydraulic conductivity (a constant), $R$ a recharge constant, and $S$ the storage coefficient, another constant. In FD solutions for the unconfined case (where the transmissivity varies with $h$) this equation is often linearized by assuming that $T$ is constant. If you ignore this linearization and continue, you get:

$$\frac{\partial}{\partial x}(Kh\frac{\partial h}{\partial x}) + R = S\frac{\partial h}{\partial t}$$

Now I don't see a problem with solving this. Couldn't one just apply the chain rule of differentiation to get...

$$\frac{\partial Kh}{\partial x}\frac{\partial h}{\partial x} +Kh\frac{\partial^2 h}{\partial x^2}+ R = S\frac{\partial h}{\partial t}$$

... and then solve each of these terms by conventional centered differences and second order central differences? I can see that this kind of solution may be problematic for implicit schemes, where you will get a squared $h$ in the first term, but why is this not done in explicit schemes?

• It can and is done in FD schemes. Some people don't need the full nonlinear equations for their modeling and so they make the linearization assumption you described. Others want the full nonlinearity and approach the problem using approaches like you describe. The thing you have to keep in mind is that solving nonlinear forms for PDEs makes stability trickier and so if you can approximate your problem with a linear version, that is more desirable. – spektr May 26 '17 at 17:03
• Thank you very much for your response. That makes sense - I can see how stability might become more tricky if you involve a nonlinear term. Then I apparently just was unfortunate enough to never stumble across a source pursuing a non-linear FD solution. Have a nice weekend! – J.Galt May 26 '17 at 17:36
• no problem! You should consider looking at how people solve the burgers equation with FD. That's pretty basic nonlinear PDE that should have some FD examples out in the internet! – spektr May 26 '17 at 18:23
• You say that your first equation originates from "Fluid Dynamics".The first term in that equation is a diffusion term which is typically not a dominant term in most fluids problems. Discretizing that term with FD leads to terms with $(\Delta x)^2$ in the denominator which in turn leads to a rather severe restriction on the time step for explicit methods. Based on my experience, that is why the equation you show is usually solved using implicit methods (e.g. Crank-Nicolson). And, as you observed, the nonlinearities complicate implicit solvers. – Bill Greene May 26 '17 at 18:41
• Assuming that T is constant also leads to a symmetric discretization of the first term. Symmetric linear systems are often much easier to solve than asymmetric ones since many methods such as conjugate gradient only work for symmetric systems. (Of course, the presence of asymmetric boundary conditions often screws this up anyway.) – nukeguy May 26 '17 at 18:56

For explicit schemes, this indeed makes sense -- there is no reason to linearize the explicit terms. However, as others have already pointed out, you get into trouble with explicit schemes for the diffusion equation because the time step must be chosen proportional to something like $$\Delta t \propto \min_x \frac{S(x) \delta x^2}{T(x)} \propto \min_x \frac{S(x) \delta x^2}{Kh},$$ i.e., you need to choose the time step very small if either $S$ becomes small somewhere in the domain, or $h$ becomes large. In addition, of course, reducing the grid size by a factor of 2 requires you to reduce the time step size by an undesirable factor of 4.
• They're unconditionally stable, but of course not unconditionally accurate: The error is still going to be something like ${\cal O}(\Delta x^p + \Delta t^q)$. So a larger time step of course means a larger error. But intuitively, for parabolic problems like yours, one would want to be able to choose $\Delta t$ larger and larger as the solution goes towards a steady state. If the increments in every time step become smaller (because you get to a steady state), you don't need to worry about accuracy as much any more. Except you can't choose $\Delta t$ large for explicit schemes. – Wolfgang Bangerth May 31 '17 at 22:03
• No, au contraire. We use implicit schemes precisely because of accuracy. Even if the system never goes into equilibrium, we still want to choose $\Delta t \propto \Delta x^{p/q}$ for accuracy reasons, not $\Delta t \propto \Delta x^2$ for stability reasons. The former almost always allows for larger time steps. – Wolfgang Bangerth Jun 3 '17 at 1:03