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.

### Newton-Raphson method

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 :

where is approximated by

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 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 and`epsilon`

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 simple`def`

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 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
- Propose a very simple test protocol for a function .
- Enhance your implementation to include backtracking inside your Newton-Raphson implementation.

### 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.

- 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 .

As a means of verification, and with the conditions above, at , one gets the following values for the total force and its Jacobian :

and

- 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 form and to apply the function.

- 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

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 :

Notice that is a vector with coordinates.

- 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`

) - 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 :

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 :

- 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.