Newton-Raphson Non-Linear Tridiagonal

jacobiannewton raphsonnonlinear systemtaylor expansion

I have the following problem where I am asked to solve a system of nonlinear equations. I am positive that I have to use Newton-Raphson with the Jacobian. My problem is that I don't fully understand the notation. This problem has $u_{i}$, $u_{i-1}$, $u_{i+1}$, are those supposed to be your $x,y,z$? If so, how do the conditions on the side of each equation apply?

The problem can be seen below:

$$\begin{align*}-3u_i+u_{i+1}&-0.5h^2u_i^3=-1.005 & i=1\\u_{i-1}-3u_i+u_{i+1}&-0.5h^2u_i^3=-0.005 &2\leq i\leq n-1\\u_{i-1}-3u_i&-0.5h^2u_i^3=-1.005 &i=n\end{align*}$$ Use $h=0.1$ and $n=30$.

The problem can also be seen here

Best Answer

I will work this out for a smaller $n = 3$, so you can see the process. Using the given iteration formula, we have

$$\begin{align}-3 u_1+u_2 -0.005 u_1^3 +1.005 = 0 \\-u_1-3 u_2+u_3 -0.005 u_2^3+0.005 = 0 \\ u_2-3 u_3 -0.005 u_3^3 +1.005 = 0 \end{align}$$

The regular Newton-Raphson method is initialized with a starting point $u_0$ and then you iterate $\tag 1 u_{n+1}=u_n-\dfrac{f(u_n)}{f'(u_n)}$

In higher dimensions, there is an exact analog. We define

$$F\left(\begin{bmatrix}u_1\\u_2\\u_3\end{bmatrix}\right) = \begin{bmatrix}f_1(u_1,u_2,u_3) \\ f_2(u_1,u_2,u_3) \\ f_3(u_1,u_2,u_3) \end{bmatrix} = \begin{bmatrix}-3 u_1+u_2 -0.005 u_1^3 +1.005 \\ -u_1-3 u_2+u_3 -0.005 u_2^3+0.005 \\ u_2-3 u_3 -0.005 u_3^3 +1.005 \end{bmatrix} = \begin{bmatrix}0\\0\\0\end{bmatrix}$$

The derivative of this system is the $3x3$ Jacobian given by

$$J(u_1,u_2,u_3) = \begin{bmatrix} \dfrac{\partial u_1}{\partial f_1} & \dfrac{\partial u_2}{\partial f_1} & \dfrac{\partial u_3}{\partial f_1}\\ \dfrac{\partial u_1}{\partial f_2} & \dfrac{\partial u_2}{\partial f_2} & \dfrac{\partial u_3}{\partial f_2} \\ \dfrac{\partial u_1}{\partial f_3} & \dfrac{\partial u_2}{\partial f_3} & \dfrac{\partial u_3}{\partial f_3}\end{bmatrix} = \begin{bmatrix} -0.015 u_1^2-3 & 1 & 0 \\ -1 & -0.015 u_2^2-3 & 1 \\ 0 & 1 & -0.015 u_3^2-3 \end{bmatrix}$$

The function $G$ is defined as

$$G(u) = u - J(u)^{-1}F(u)$$

and the functional Newton-Raphson method for nonlinear systems is given by the iteration procedure that evolves from selecting an initial $u^{(0)}$ and generating for $k \ge 1$ (compare this to $(1)$),

$$u^{(k)} = G(u^{(k-1)}) = u^{(k-1)} - \dfrac{F(u^{(k-1)})}{J(u^{(k-1)})} = u^{(k-1)} - J(u^{(k-1)})^{-1}F(u^{(k-1)}).$$

We can write this as

$$\begin{bmatrix}u_1^{(k)}\\u_2^{(k)}\\u_3^{(k)}\end{bmatrix} = \begin{bmatrix}u_1^{(k-1)}\\u_2^{(k-1)}\\u_3^{(k-1)}\end{bmatrix} + \begin{bmatrix}v_1^{(k-1)}\\v_2^{(k-1)}\\v_3^{(k-1)}\end{bmatrix}$$

where

$$\begin{bmatrix}v_1^{(k-1)}\\v_2^{(k-1)}\\v_3^{(k-1)}\end{bmatrix}= -\left(J\left(u_1^{(k-1)},u_2^{(k-1)}, u_3^{(k-1)}\right)\right)^{-1}F\left(u_1^{(k-1)},u_2^{(k-1)},u_3^{(k-1)}\right)$$

Using the starting vector

$$u^{(0)} = \begin{bmatrix}u_1^{(0)}\\u_2^{(0)}\\u_3^{(0)}\end{bmatrix} = \begin{bmatrix}0\\0\\0\end{bmatrix}$$

You end up with (the iteration converges very quickly in this case)

$$\begin{bmatrix}u_1\\u_2\\u_3\end{bmatrix} = \begin{bmatrix} 0.33549261976670497\\ 0.00166666665895062\\ 0.33549261976670497 \end{bmatrix}$$

Of course, for a different starting vector you may get a different solution or no solution at all.

Note: It is worth noting that the way you code this up should provide a solution for any sized $n$, that is, your code should easily generalize regardless of $n$, while adhering to the constraints of the numerical methods you are using, and that is likely the point of this exercise. This means that your code should be able to define everything you need from the iteration formulas themselves. You just enter $n$ and the rest should be automated.

Update for the case $n = 4$, the iteration formula gives

$$\begin{align}-3 u_1+u_2 -0.005 u_1^3 +1.005 = 0 \\ u_1 - 3 u_2 + u_3 -0.005 u_2^3+0.005 = 0 \\ u_2-3 u_3 + u_4 -0.005 u_3^3 +0.005 = 0 \\ u_3 - 3 u_4 -0.005 u_4^3 + 1.005 = 0\end{align}$$

Related Solutions

[Math] Multivariate Newton-Raphson in R language (equations that contains integral)

pp.13-16 here discuss a library function that does what you need to use Newton-Raphson, the multiroot function in the rootSolve package. The compulsory arguments of multiroot are a function f, which for your purposes will send a 2D vector to a 2D vector, and an initial value for its argument so you can begin the iteration. The real challenge is creating the function f, whose first and second components should be the values of F_a and F_b. So you want something like:

F_a <- function(a, b) integrate(function(z) some_function_f(a, b), -Inf, Inf)$value-a

F_b <- function(a, b) integrate(function(z) some_function_g(a, b), -Inf, Inf)$value-b

f <- function(v)
{
    a <- v[1]
    b <- v[2]
    c(F_a(a, b), F_b(a, b))
}

solve_for_a_and_b <- function(original_a, original_b) rootSolve::multiroot(f, c(original_a, original_b))

Jacobian of a system of equations including an ODE for Newton-Raphson

The general idea of Newton-like methods is that you approximate the equation by some linearization. This linearization does not need to be very exact to still get a reasonably fast convergence. Here just iterating the stage equations as a fixed-point iteration gives linear convergence with a factor $O(h)$. If the linearization is also $O(h)$ exact, then the resulting Newton-like iteration will still be linear, but now with a factor $O(h^2)$. Depending on how the initial guess, the predictor, is computed, this can result in a drastic reduction of required corrector steps in the Newton-like iteration to get to an error level that is below the method step truncation error.

$\newcommand{\D}{\mathit{\Delta}}$ So if $J$ is the Jacobian at $y_n$, then you can decompose $f(y+\D y)=f(y)+J\D y+R(\D y)$. For the given 2-stage method this gives the system $$ \vec k_1-hJ(B_{11}\vec k_1+B_{12}\vec k_2)=f(\vec y_n)+R(h(B_{11}\vec k_1+B_{12}\vec k_2)) \\ \vec k_2-hJ(B_{21}\vec k_1+B_{22}\vec k_2)=f(\vec y_n)+R(h(B_{21}\vec k_1+B_{22}\vec k_2)) $$ In an iteration step the left side is a linear system for the new values, while the right side is computed from the old values. This form demonstrates that the right side is constants plus small superlinear corrections, in the practical computation the right side will be computed as $f(y)+R(\D y^{old})=f(y+\D y^{old})-J\D y^{old}$. \begin{align} \pmatrix{I_2-hB_{11}J&-hB_{12}J\\-hB_{21}J&I_2-hB_{22}J} \pmatrix{\vec k_1\\\vec k_2} = \pmatrix{f(y+\D_1 y^{old})-J\D_1 y^{old}\\f(y+\D_2 y^{old})-J\D_2 y^{old}} \end{align} Using a Kronecker product, the matrix on the left could also be written as $I_4-hB\otimes J$.

Best Answer

Related Solutions

[Math] Multivariate Newton-Raphson in R language (equations that contains integral)

Jacobian of a system of equations including an ODE for Newton-Raphson

Related Question