Efficient solution of large systems of non linear algebraic equations

Question

I have quite simple problem of FEM solution of 1D differential non-linear equation. Altough the problem itself is simple, the numerical solution of the arising system of non-linear algebraic equations takes hours. The system is tridiagonal type, typical for 1D FEM problems. At the moment I am using Picard or Newton methods. And those are the only methods I have experience with. However they do not seem to give good results. Are there other methods, modifications to those or solvers anyone could recommend?

EDIT: The problem is solution of 1D diffusve wave equation. The details can be found here. The system is

$$\textbf{M}\cdot{\textbf{x}}={\textbf{r}}$$ with $$M_{i,i-1}=(1-\omega)\cdot \frac{\Delta x}{2}+\theta\cdot \Delta t \cdot \frac{1}{2}\left( \frac{1}{n}\frac{(x_{j-1}-z_{j-1})^{5/3}}{ \sqrt{\frac{x_{j}-x_{j-1}}{\Delta x}}}+ \frac{1}{n}\frac{(x_{j}-z_{j})^{5/3}}{ \sqrt{\frac{x_{j}-x_{j-1}}{\Delta x}}}\right)$$

$$M_{i,i}=\omega\cdot \Delta x+\theta\cdot \Delta t \cdot \frac{1}{\Delta x}\left [ \frac{1}{2}\left( \frac{1}{n}\frac{(x_{j-1}-z_{j-1})^{5/3}}{ \sqrt{\frac{x_{j}-x_{j-1}}{\Delta x}}}+ \frac{1}{n}\frac{{}(x_{j}-z_{j})^{5/3}}{ \sqrt{\frac{x_{j}-x_{j-1}}{\Delta x}}}\right) +\frac{1}{2}\left( \frac{1}{n}\frac{(x_j-z_j)^{5/3}}{ \sqrt{\frac{x_{j+1}-x_j}{\Delta x}}}+ \frac{1}{n}\frac{{}(x_{j+1}-z_{j+1})^{5/3}}{ \sqrt{\frac{x_{j+1}-x_j}{\Delta x}}}\right) \right] $$

$$M_{i,i+1}=(1-\omega)\cdot \frac{\Delta x}{2}+\theta\cdot \Delta t \cdot \frac{1}{2}\left( \frac{1}{n}\frac{(x_j-z_j)^{5/3}}{ \sqrt{\frac{x_{j+1}-x_j}{\Delta x}}}+ \frac{1}{n}\frac{{}(x_{j+1}-z_{j+1})^{5/3}}{ \sqrt{\frac{x_{j+1}-x_j}{\Delta x}}}\right)$$

With $\textbf{r}$ constant. So concluding the non-linearity is quite strong as there is 1 by square root and additionally unknown value to the power of $5/3$. For solution of system of linear equations I use lu decomposition for sparse matrices. I do not consider those systems large, this word was added by moderator. The systems are of size 1000x1000.

Answer 1

Looking to the paper you referred to in EDIT, my first guess for possible problems in your nonlinear solver is that in fact the nonlinear system you present is different (simplified) than the one given in the paper.

In your description you do not mention what to do if during iterations one obtains e.g. $x_{j+1}-x_j \approx 0$ when you divide by zero or if $x_j < z_j$ when you get complex numbers. Moreover following the paper you should have instead of $\sqrt{\frac{x_{j+1}-x_j}{\Delta x}}$ the term $\sqrt{\frac{|x_{j+1}-x_j|}{\Delta x}}$. The nonlinearity in your algebraic system comes from the approximation of nonlinear diffusion coefficient in PDE where the cases I mentioned are in the paper defined differently (``regularized'').

If you do follow the paper also in these definitions (it means you simply did not describe it in your question), you should really specify more precisely what you do not like in the behavior of your solvers, but such ``piecewise'' defined non-linearity as mentioned in the paper can deteriorate the behavior of nonlinear iterations.

Nevertheless, I would expect that if you solve a problem where those cases can not happen e.g. a propagation of a wave given initially by $h(x,0)=2+sin(x)$ for $z \equiv 0$ and boundary conditions $h(-\pi/2,t)=1$ and $h(\pi/2,t)=3$, you really should get a quadratic convergence of Newton solver for enough small time step $\Delta t$. I use the same notation as in the referred paper, it means the system is in fact $M h = r$ and $x$ is the space variable, so the discrete unknowns are $h_j$ and not $x_j$.