Finite difference equations representing semilinear elliptic PDE

Question

I recently asked a question pertaining to the appliciation of Jacobi's method to a semilinear elliptic PDE (Poisson's equation)

$$ \nabla^2u = -\rho~e^{-u} $$

A more efficient method like the Bi conjugate gradient stabilised method was recommended. I have tested this method out and it is indeed much faster. But I am unsure of what the matrix representation of a semilinear system would look like. For an ordinary linear PDE like

$$ \nabla^2u=-\rho $$

it looks like $$ \frac{1}{h^2}\left( \begin{array}{ccc} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} \rho_1+g \\ \rho_2 \\ \rho_3 \\ \rho_4 \\ \rho_5+g\end{array} \right) $$ where $g$ is the Dirichlet boundary condition.

My question: What would the corresponding matrix representation of the set of simultaneous equations for the semilinear case look like? I'm guessing something like

$$ \frac{1}{h^2}\left( \begin{array}{ccc} 2 & -1 & & & \\ -1 & 2 & -1 & & \\ & -1 & 2 & -1 & \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \\ \end{array} \right)\left( \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5\end{array} \right) = \left( \begin{array}{c} \rho_1e^{-u_1}+g \\ \rho_2e^{-u_2} \\ \rho_3e^{-u_3} \\ \rho_4e^{-u_4} \\ \rho_5e^{-u_5}+g\end{array} \right) $$

But this doesn't leave me with all $u$ values in a single vector.

Would it make sense to do something like:

1) Solve the linear case $$ \nabla^2u = -\rho $$ 2) Use the resultant $u$ to construct a new linear case $$ \nabla^2u_i = -\rho~C $$ where $$ C = e^{-u_{old}} $$ 3) Repeat step 2 until self-consistency is reached.

Answer 1

Since you don't end up with a linear system, you have to use a nonlinear root finder. The most common choice is Newton's method, because it converges quite quickly (if your initial guess is good enough). You've probably seen Newton's method for a single variable in a single equation. There is an analogous method for systems of equations.

Given a vector function $f(\mathbf{u}^*) = 0$, and an initial guess $\mathbf{u}^0$, we can approximate the root $\mathbf{u}^*$ by solving the sequence of linear systems $$ J(\mathbf{u}^k)\Delta\mathbf{u}^k = -f(\mathbf{u}^k),$$ where $\Delta \mathbf{u}^k = \mathbf{u}^{k+1} - \mathbf{u}^{k}$ and $J(\mathbf{u}^k)$, called the Jacobian, is defined as

$$ J(\mathbf{u}^k) = \begin{pmatrix} \frac{\partial f_1}{\partial u_1}(\mathbf{u}^k) & \cdots & \frac{\partial f_1}{\partial u_n}(\mathbf{u}^k)\\ \vdots & \ddots & \vdots\\ \frac{\partial f_n}{\partial u_1}(\mathbf{u}^k) & \cdots &\frac{\partial f_n}{\partial u_n}(\mathbf{u}^k)\end{pmatrix}.$$

In your case $f(\mathbf{u})$ is given by

$$ f(\mathbf{u}) = \frac{1}{h^2}A\mathbf{u} - \mathbf{b}(\mathbf{u})$$ where $A$ is your finite difference matrix and $\mathbf{b}(\mathbf{u})$ is your right hand side. The Jacobian is given by

$$ J(\mathbf{u}) = \frac{1}{h^2}A + \begin{pmatrix} \rho_1 e^{-u_1}\\ \mathbf{0} & \ddots & \mathbf{0} \\ & & \rho_n e^{-u_n}\end{pmatrix}.$$

Your initial guess could come from solving the homogeneous system for example. Then given $\mathbf{u}^0$, we can evaluate $J(\mathbf{u}^0)$ and $f(\mathbf{u}^0)$ and solve for $\Delta \mathbf{u}^0$ using $$ J(\mathbf{u}^0)\Delta\mathbf{u}^0 = -f(\mathbf{u}^0).$$

We then find $\mathbf{u}^1$ from the definition of $\Delta \mathbf{u}^0$. We can repeat the process to find $\mathbf{u}^2, \mathbf{u}^3, \cdots$, stopping once $|\Delta \mathbf{u}^k|$ is below some tolerance.

Answer 2

@H.H. already gives a good answer on how to implement things for the Newton iteration. If you want to get an overview of different approaches to address nonlinear problems, take a look at lectures 31.5 and following here: http://www.math.colostate.edu/~bangerth/videos.html .