The goal of this project is to program algorithms dedicated to the research of roots of systems of non-linear equations. The method promoted here is the Newton-Raphson algorithm, and the goal is to evaluate the assets and liabilities of such a solution. This shall be done by testing the method in different settings. In the following, you will be proposed a list of applications where this algorithm is necessary. You must program at least 2 of these applications, and then write a summary of your experiments as a conclusion.
- Chapter 6 of « Methodes Numeriques, Algorithmes, analyse et applications ». The book is available on my enseirb account in ~mfaverge/public_html/IS104/Methodes-Numeriques-Algorithmes-analyse-et-applications.pdf
- The slides made to help presenting this chapter here. They are also available on my account at enseirb alongside with the book.
Newton-Raphson method


![f : \left[\begin{array}{c} x_1\\\vdots\\x_n\end{array}\right] \rightarrow \left[\begin{array}{c} f_1(x_1,\dots,x_n)\\\vdots\\ f_m(x_1,\dots,x_n)\end{array}\right]](wp-content/cours/IS104/images/f:left[begin{array}{c}x_1vdotsx_nend{array}right]rightarrowleft[begin{array}{c}f_1(x_1dotsx_n)vdotsf_m(x_1dotsx_n)end{array}right].png)
The function is called the Jacobian matrix of
. The idea of the Newton-Raphson algorithm consists in considering that, given a point
, the best direction leading to a root of the function is the direction given by the slope of
at
. In following this direction, the algorithm assumes that it moves closer to a root of
.
Given a vector representing the current position, the algorithm computes
, a vector representing the value of
at
, as well as
the Jacobian matrix of the derivatives of
at
. Now, the next position
+
is chosen such that :
=0.png)
.png)
+underbrace{H(U)}_{textrm{matrix}}timesunderbrace{V}_{textrm{vector}}.png)
Then, is replaced by
, and this operation is repeated until either
converges or a maximal number of iterations is reached.
-
- Write the (very simple) equation linking
,
and
.
- Write the (very simple) equation linking
-
- Write the Newton-Raphson method in a generic way :
function U = Newton_Raphson(f, J, U0, N, epsilon)
… where
f
is the function under study,J
is its Jacobian matrix,U0
is the starting position of the algorithm,N
is the maximal number of steps allowed during the algorithm andepsilon
measures the convergence of the algorithm.Remarks :- One requirement is to be able to compute the values taken by the function
and all its derivatives in a priori any point of the plane. Therefore, it is necessary to specify
and its derivatives in the form of functions. Python allows programming using functional parameters. These parameters may be provided by functions defined either with the keyword
lambda
, or with simpledef
definitions. - The algorithm requires the resolution of several linear systems with matrices
that are possibly singular (i.e
. Prefer the function
numpy.linalg.lstsq
in order to avoid such numerical problems. - One of the drawbacks of this algorithm is its tendency to diverge in numerous cases. It is therefore imperative to limit the number of steps taken, and to detect to what extent the algorithm has converged.
- One requirement is to be able to compute the values taken by the function
- Propose a very simple test protocol for a function
.
- Enhance your implementation to include backtracking inside your Newton-Raphson implementation.
- Write the Newton-Raphson method in a generic way :
Computation of the Lagrangian points
-
- Write the code to represent the following kinds of forces :
-
- an elastic force (with zero natural length)
- a centrifugal force
- a gravitational force
Remarks :- Is is strongly recommended to use functional programming techniques to represent the forces and their Jacobian matrices.
- For each of these cases, it should be possible to parameterize the force by a constant
symbolizing its intensity, as well as a central point
from which the force is issued.
-
- Write the code to represent the following kinds of forces :
-
- Use the Newton-Raphson method to obtain the equilibrium points in the following case :
-
- Two gravitational forces with respective coefficients
(resp.
) originating from
(resp.
);
- And a centrifugal force centered on the barycenter of the two masses, with coefficient
.
- Two gravitational forces with respective coefficients
As a means of verification, and with the conditions above, at
, one gets the following values for the total force and its Jacobian :
and
-
- Use the Newton-Raphson method to obtain the equilibrium points in the following case :
-
- Is it possible to obtain the Lagrangian points ?
Bairstow method for polynomial factorization



-
- Describe (literally) the function whose zero is computed in dimension 2. In particular, define its domain of definition and its range, and explain the role of the variables
,
,
and
. What are the zeroes of this function ?
- Write the code of this function by reusing the function
numpy.polydiv
(which computes simultaneously the quotient and the remainder of the division). In order to test this function, it suffices to construct the adequate polynomial of the formand to apply the function.
- Describe (literally) the function whose zero is computed in dimension 2. In particular, define its domain of definition and its range, and explain the role of the variables
-
- In the Numerical Recipes, find the 4 partial derivatives of the previous function. Write the equations defining these derivatives, and program them.
-
- Write the Bairstow algorithm, and test it using a method of your choice.
Electrostatic equilibrium
![[-1, 1]](wp-content/cours/IS104/images/[-11].png)


![[-1, 1]](wp-content/cours/IS104/images/[-11].png)
=sum_{i=1}^Nlog|x_i+1|+log|x_i-1|+frac{1}{2}sum_{j=1jnei}^Nlog|x_i-x_j|.png)
The equilibrium positions are found by minimizing or maximizing this energy. In order to do this, it is necessary to solve a non-linear system of equations that is equal to :
![\nabla E(x_1, x_2, \cdots , x_N) = \displaystyle\left[ \begin{array}{c} \frac{\partial E(x_1,\dots,x_N)} {\partial x_i} \\ \end{array}\right] = 0](wp-content/cours/IS104/images/nablaE(x_1x_2cdotsx_N)=displaystyleleft[begin{array}{c}frac{partialE(x_1dotsx_N)}{partialx_i}end{array}right]=0.png)
Notice that
.png)

-
- Compute the Jacobian of
.
- Compute the Jacobian of
-
- Use Newton’s method in order to solve this equation. Plot the points and the real axis.
Do these solutions resemble to the roots of the derivative of the Legendre polynomials ? (cf.numpy.polynomial.legendre
)
- Use Newton’s method in order to solve this equation. Plot the points and the real axis.
-
- Does the solution correspond to a maximum or a minimum of the energy ?
The wave equation

where denotes the derivative of
with respect to
.
is the « scaled » elevation of the surface of the wave and is a function of the time
and the space
. The true elevation is equal to
, the parameter
will be chosen less than 1 (e.g. 0.1). We will reduce the computations on the interval
and impose periodic boundary conditions.
Let be a vector of size
, discretizing
on a regular subdivision
of the interval
(
). We are going to compute a sequence of vector
representing the temporal evolution of the wave. The initial condition will be set to a gaussian function :
with
. As for the heat equation, a finite difference scheme is used to solve this equation. In the following, we describe an algorithm to compute
as a function of
.
Periodic boundary conditions induce the following conditions :

We denote the application which transforms the vector
into the vector :
![[ G(U) ]_i = \frac{U_{i+1} - U_{i-1}}{2 \Delta x} + \frac{\varepsilon}{4 \Delta x} (U_{i+1} - U_{i-1})(U_{i-1} + U_i + U_{i+1}) + \frac{\varepsilon}{12 \Delta x^3} (U_{i+2} - 2 U_{i+1} + 2 U_{i-1} - U_{i-2})](wp-content/cours/IS104/images/[G(U)]_i=frac{U_{i+1}-U_{i-1}}{2Deltax}+frac{varepsilon}{4Deltax}(U_{i+1}-U_{i-1})(U_{i-1}+U_i+U_{i+1})+frac{varepsilon}{12Deltax^3}(U_{i+2}-2U_{i+1}+2U_{i-1}-U_{i-2}).png)
This application is non-linear and can not be represented as a matrix. We have the following evolution system : . In order to solve this evolution equation, we consider the following equation :
=0.png)
-
- Compute and draw
.
- Write the function whose roots you must find. Compute the Jacobian matrix of this function.
- Compute
as a function of
with the Newton-Raphson method. Draw the results. As an example, it should be possible to generate a movie like this one.
- Compute and draw