Page MenuHomec4science
Files F60948512

NonLinearEquationsLib.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 13:36

NonLinearEquationsLib.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 30
@author: Simone Deparis
"""
import numpy as np
def FixedPoint( phi, x0, a, b, tol, nmax ) :
'''
FixedPoint Find the fixed point of a function by iterative iterations
FixedPoint( PHI,X0,a,b, TOL,NMAX) tries to find the fixedpoint X of the a
continuous function PHI nearest to X0 using
the fixed point iterations method.
[a,b] : if the iterations exit the interval, the method stops
If the search fails an error message is displayed.
Outputs : [x, r, n, x_sequence]
x : the approximated fixed point of the function
r : the absolute value of the residual in X : |phi(x) - x|
n : the number of iterations N required for computing X and
x_sequence : the sequence computed by Newton
'''
# Initial values
n = 0
xk = x0
# initialisation of loop components
# increments (in abs value) at each iteration
inc = []
# in case we wish to plot the sequence
x = [x0]
# diff : last increment,
diff = tol + 1 # initially set larger then tolerance
# Loop until tolerance is reached
while ( diff >= tol and n <= nmax and (xk >= a) and (xk <= b) ) :
# call phi
xk1 = phi(xk)
# increments
diff = np.abs(xk1 - xk)
# prepare the next loop
n = n + 1
xk = xk1
x.append(xk)
# Final residual
rk1 = np.abs(xk1 - phi(xk1))
# Warning if not converged
if n > nmax :
print('FixexPoint stopped without converging to the desired tolerance ')
print('because the maximum number of iterations was reached')
return xk1, rk1, n, np.array(x)
def plotPhi (a,b,phi,label='$\phi$', N = 100) :
'''
simple plot of fonction with bisectrice of first quadrant
useful for preparing grpahics of fixed point iterations
[a,b] : interval, used for both x and y axis
phi : funciton to plot
label : label for the fonction, usually '$\phi$'
N : number of points for plotting
'''
import matplotlib.pyplot as plt
z = np.linspace(a,b,N)
plt.plot(z,phi(z),'k-')
plt.plot([a,b],[a,b],':', linewidth=0.5)
plt.xlabel('x'); plt.ylabel(label);
# Plot the x,y-axis
plt.plot([a,b], [0,0], 'k-',linewidth=0.1)
plt.plot([0,0], [a,b], 'k-',linewidth=0.1)
plt.legend([label, 'y=x'])
plt.title('Graph de la fonction ' + label)
def plotPhiIterations (x) :
# plot the graphical interpretation of the Fixed Point method
import matplotlib.pyplot as plt
plt.plot([x[0],x[0]], [0,x[1]], 'g:')
for k in range(x.size-2) :
plt.plot([x[k],x[k+1]], [x[k+1],x[k+1]], 'g:')
plt.plot([x[k+1],x[k+1]], [x[k+1],x[k+2]], 'g:')
k=x.size-2
plt.plot([x[k],x[k+1]], [x[k+1],x[k+1]], 'g:')
# Putting a sign at the initial point
deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0]
if (x[1] < 0) :
plt.annotate("$x_0$", (x[0], 0.02*deltaAxis) )
else :
plt.annotate("$x_0$", (x[0], -0.05*deltaAxis) )
def plotNewtonIterations (a,b,f,x,N=200) :
# plot the graphical interpretation of the Newton method
import matplotlib.pyplot as plt
z = np.linspace(a,b,N)
plt.plot(z,f(z), 'b-', x,f(x), 'rx')
# Putting a sign at the initial point
deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0]
# plot the graphical interpretation of the Newton method
plt.plot([x[0],x[0]], [0,f(x[0])], 'g:')
for k in range(x.size-1) :
plt.plot([x[k],x[k+1]], [f(x[k]),0], 'g-')
plt.plot([x[k+1],x[k+1]], [0,f(x[k+1])], 'g:')
# Putting a sign at the initial point
if (f(x[k]) < 0) :
plt.annotate("$x_"+str(k)+"$", (x[k], 0.02*deltaAxis) )
else :
plt.annotate("$x_"+str(k)+"$", (x[k], -0.05*deltaAxis) )
plt.ylabel('$f$'); plt.xlabel('$x$');
# Plot the x,y-axis
plt.plot([a,b], [0,0], 'k-',linewidth=0.1)
plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)
plt.legend(['$f$','($x_k$,$f(x_k)$)'])
def plotBisectionIterations (a,b,f,x,N=200) :
# plot the graphical interpretation of the Bisection method
import matplotlib.pyplot as plt
def putSign(y,f,text) :
if (f(y) < 0) :
plt.annotate(text, (y, 0.02*deltaAxis) )
else :
plt.annotate(text, (y, -0.05*deltaAxis) )
z = np.linspace(a,b,N)
plt.plot(z,f(z), 'b-')
plt.plot([a,a], [0,f(a)], 'g:')
plt.plot([b,b], [0,f(b)], 'g:')
# For putting a sign at the initial point
deltaAxis = plt.gca().axes.get_ylim()[1] - plt.gca().axes.get_ylim()[0]
putSign(a,f,"a")
putSign(b,f,"b")
for k in range(x.size) :
plt.plot([x[k],x[k]], [0, f(x[k])], 'g:')
putSign(x[k],f,"$x_"+str(k)+"$")
# Plot the x,y-axis
plt.plot([a,b], [0,0], 'k-',linewidth=0.1)
plt.plot([0,0], [np.min(f(z)),np.max(f(z))], 'k-',linewidth=0.1)
plt.ylabel('$f$'); plt.xlabel('$x$');
return

Event Timeline

c4science · Help