For the numerical solution of Reynolds equations (a non-linear partial differential equation), the Newton-Raphson method is generally proposed.

After getting algebraic equations from a finite difference discretization, the Newton-Raphson method is applied to those non-linear algebraic equations. This yields a square system of linear equations.

In the case of a two-dimensional mesh with $m\times n$ grid points ($m\neq n$) system, is it necessary to have $n>m$?

Also, how do we solve those linear equations?

  • $\begingroup$ The plural of "matrix" is "matrices". $\endgroup$
    – joriki
    May 4 '12 at 11:57
  • 1
    $\begingroup$ Since it's pretty obvious answers/comments will no longer help the original author, I am wondering if this question is worth salvaging. The wording of it makes figuring out what is actually being asked fairly hard to discern. I am assuming he/she is referring to a structured rectangular mesh, and the wording regarding newton-raphson "linearizing" the non-linear equations is very confusing. Then the question in the final line appears to ask about methods for matrix factorization/solutions. $\endgroup$ Sep 27 '12 at 2:08

For solving nonlinear partial equations, say $-\Delta y + f(y) = 0$, there are two common approaches:

  • Fixed point iteration: Pick $y^0$ and for $k=1,\dots$ solve $$ -\Delta y^{k+1} + f(y^k) = 0.$$ (That is, you replace $y$ by $y^{k+1}$ in all linear terms and by $y^k$ in all nonlinear terms.)

  • Newton's method: Pick $y^0$ and for $k=1,\dots$ solve $$ -\Delta \delta y + f'(y^k)\delta y = -(\Delta y^k+f(y^k)) \tag{1} $$ and set $y^{k+1} = y^k + \delta y$.

Which to use depends on the nonlinear term $f$; in the context of Reynolds equation, $f$ is quadratic, and this means that Newton's method is a good choice.

Now to the main question: How to actually compute $\delta y$. The easiest (and, at least for me, least error-prone way) is to apply your finite-difference discretization to $(1)$ and solve the discretized problem like you would any linear PDE.

Alternatively, let's say you already have a discretization of the nonlinear problem and want to solve the resulting nonlinear finite-dimensional system. From your question, I take it that you have a two-dimensional domain which you have discretized using $N:=mn$ grid points, resulting in matrices $K_h\in\mathbb{R}^{N\times N}$ (representing the linear part) and $F_h(y_h)\in\mathbb{R}^{N\times N}$ (representing the nonlinear part), where $y_h\in\mathbb{R}^N$ is the discrete solution. Since the dimensions only depend on the product $N$, this means that Newton's method does not care whether $n>m$ or not (your discretization will, of course). Applying Newton's method to $K_h y_h + F_h(y_h) = 0$, you have to solve for $\delta y_h \in \mathbb{R}^N$ in $$(K_h + \nabla F_h(y_h^k)) \delta y_h = - K_h y_h^k - F_h(y_h^k),$$ where $\nabla F_h(y_h)\in\mathbb{R}^{N\times N}$ is the Jacobi matrix of $F_h$ with entries $(\nabla F_h(y))_{ij} = \frac{\partial F_h(y_h)_i}{\partial (y_h)_j}$. This is a standard linear system, which you can solve with any of the usual methods. There is one additional point: If you pick an iterative method, you can exploit the fact that for $y^k$ far from the solution $y$, your linearization error will be large and likely dominate the approximation error, and thus there is no need to solve the linear system very accurately in the beginning. (Of course, you then need to increase the accuracy as you get closer to $y$.) This is known as an inexact Newton(-Krylov) method.

  • $\begingroup$ Just noticed a few typos: the right-hand side of (1) should be $-f(y^k) + \Delta y^k$ and $F_h(y_h) \in \mathbb{R}^N$ is not a matrix, of course. (I don't want to edit the post and bump the question.) $\endgroup$ Aug 18 '14 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy