Commit 2ab2a0af authored by Swicech's avatar Swicech
Browse files

Update Newton_s_method_for_systems_of_equations_assignment__with_solutions_.ipynb

parent 82894b65
%% Cell type:markdown id: tags:
# Newton's method for systems of equations
Newton's method is used to solve a *general* system of equations, written in vector form
$$\mathbf{f}(\mathbf{x}) = \mathbf{0}$$
Given an initial guess $\mathbf{x}_{0}$, Newton's method is used to compute successively better approximations to the solution of the above equation by the following iterative formula
$$\mathbf{x}_{k+1} =\mathbf{x}_{k} + \mathbf{f}'(\mathbf{x_k})^{-1}\mathbf{f}(\mathbf{x_k}),$$
$$\mathbf{x}_{k+1} =\mathbf{x}_{k} - \mathbf{f}'(\mathbf{x_k})^{-1}\mathbf{f}(\mathbf{x_k}),$$
where the Jacobain derivative is
$$\mathbf{f}'(\mathbf{x}) =
\left(\begin{array}{cc}
\frac{\partial f_{0}}{\partial x_0} & \frac{\partial f_{0}}{\partial x_1} \\
\frac{\partial f_{1}}{\partial x_0} & \frac{\partial f_{1}}{\partial x_1} \\
\end{array}\right).$$
Your stopping condition should be $\|\mathbf{f}(\mathbf(x_k))\|<tol$ or $\|\mathbf{x}_k-\mathbf{x}_{k+1}\|<tol$ (or both). Where $\|\mathbf{x}\|$ is the 2-norm of the vector $\mathbf{x}$ (i.e., Pythagoras' theorem).
Here we have written the vectors with 0-based index to make it consistent with how you would index them in Python.
### Q1)
In this assignment, we will use Newton's method for solving the following set of algebraic equations:
\begin{gather}
x_1−x^3_0=0 \\
x^2_0+x^2_1−1=0
\end{gather}
#### 1a)
Write a function called **f(x)**, that accepts a vector $\mathbf{x} = (x_0,x_1)^T$ and returns a vector that corresponds to $\mathbf{f}(\mathbf{x})$ such that the above set of equations can be written in the form $\mathbf{f}(\mathbf{x}) = \mathbf{0}$.
If you have written your code correctly, then
import numpy as np
x = np.array([1,1])
print(f(x))
should return
[0 1]
%% Cell type:code id: tags:
``` python
```
%% Output
[0 1]
%% Cell type:markdown id: tags:
#### 1b)
Write a function called df(x) that accepts a vector $\mathbf{x} = (x_1,x_2)^T$ and returns the jacobian matrix for your above definition of f(x). Note: you have to do the derivatives by hand! (you send in the vector to the function and compute the calculated derivatives inside the function.)
If you have written your code correctly, then
x = np.array([1,1])
print(df(x))
should return
[[-3 1]
[ 2 2]]
%% Cell type:code id: tags:
``` python
```
%% Output
[[-3 1]
[ 2 2]]
%% Cell type:markdown id: tags:
#### 1c)
To implement Newton's method,
$$\mathbf{x}_{k+1} =\mathbf{x}_{k} + \underbrace{\mathbf{f}'(\mathbf{x_k})^{-1}\mathbf{f}(\mathbf{x_k})}_{\mathbf{v}}$$
we have to solve for the vector
$$\mathbf{v} = \mathbf{f}'(\mathbf{x_k})^{-1}\mathbf{f}(\mathbf{x_k})$$
To do this, instead of inverting the matrix $\mathbf{f}'(\mathbf{x_k})^{-1}$, we will directly solve for $\mathbf{v}$ by solving the following linear system
$$\mathbf{f}'(\mathbf{x_k})\mathbf{v} = \mathbf{f}(\mathbf{x_k}).$$
Write a function called vsolve(x,df,f), that accepts a vector x, the Jacobian df and the function f and returns a vector v that solves the above linear system.
Note: you can use a built-in function for solving linear systems such as scipy.linalg.solve
If you have written your code correctly, then
x = np.array([1,1])
print(vsolve(x,df,f))
should return
[0.125 0.375]
%% Cell type:code id: tags:
``` python
```
%% Output
[0.125 0.375]
%% Cell type:markdown id: tags:
#### 1d)
Write a function called **newton(x,f,df,tol)**, that accepts a vector **x**, a vector-valued function **f**, it's Jacobian derivative **df**, a tolerance, and returns a solution **x**, that solves **f(x)=0**.
Run your code with the initial guess $\mathbf{x} = (0,0)^T$ and $tol = 10^{-15}$.
If your code has been written correctly, then it should return a value $\mathbf{x}$ that solves $\mathbf{f}(\mathbf{x})$ to within your tolerance.
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
# test
tol = 1e-15
x0 = np.array([0, 0])
x = newton(x,f,df,tol)
print(la.norm(f(x))<tol)
```
%% Output
True
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment