Solve nonlinear reaction-diffusion PDE using implicit finite difference method

Question

I want to solve a reaction-diffusion PDE, for example
$$
\frac{\partial s}{\partial t} = D\frac{\partial^2s}{\partial x^2}+Ks^2+f(x), x, t\in[0,1]\\
s(0,t)=s(1,t)=s(x,0)=0
$$
where $s$ is a function of both $x$ and $t$. First I tried explicit finite difference method but this does not ensure stability and causes a problem. Therefore, I want to apply implicit method. For example, backward Euler method applied to a uniform grid gives
$$
\frac{s^{n+1}_i-s^n_i}{\Delta t} = D\frac{s^{n+1}_{i+1}-2s^{n+1}_i+s^{n+1}_{i-1}}{(\Delta x)^2} + K(s^{n+1}_i)^2 + f(x_i)
$$
where superscript denotes time $t=t_n$, subscript denotes space $x=x_i$.
If there is no nonlinear term $(s^{n+1}_i)^2$, I can solve this by matrix inversion. However, this nonlinear term cause trouble. Also, there are two unknowns $s_{i+1}^{n+1}$ and $s_i^{n+1}$, so Newton's method is not applicable. What is the best way to solve this?

Answer 1

One typical idea is to apply Newton's method or a Newton-like method. Essentially, you put all terms that are linear in the next state $s^{n+1}$ on the left side, and add some additional term to (locally) compensate for the non-linear terms. $$ s^{n+1}_i-D\frac{Δt}{(Δx)^2}(s^{n+1}_{i+1}-2s^{n+1}_i+s^{n+1}_{i-1})-2KΔt(a_is^{n+1}_i)=s^n_i+Δt\Bigl(K(s^{n+1}_i-a_i)^2 - Ka_i^2+ f(x_i)\Bigr) $$ Applying the inverse linear matrix of the left side to this system results in a fixed-point equation $s^{n+1}=G(s^{n+1})$. One easy choice is to fill the $a_i$ with the previous values of $s^{n+1}_i$, which removes the first term of the right side. This would give the usual Newton method where the linear system changes in each step. One could also keep the $a_i$, and thus the matrix and its factorization, constant. The $a_i$ would be filled with the initial value for the $s^{n+1}_i$, obtained via Euler step or extrapolation, or even set to $a_i=s^n_i$. This will all converge with different speeds, slower in each variant, but each still faster than the naive fixed-point iteration $$ s^{n+1}_i=G_i(s^{n+1})=s^n_i+D\frac{Δt}{(Δx)^2}(s^{n+1}_{i+1}-2s^{n+1}_i+s^{n+1}_{i-1})+Δt\Bigl(K(s^{n+1}_i)^2 + f(x_i)\Bigr). $$