Page MenuHomec4science
Files F60637032

NonLinearEquationsLib.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 1, 15:26

NonLinearEquationsLib.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan 30
@author: Simone Deparis
"""
import numpy as np
def Bisection(a,b,fun,tolerance,maxIterations) :
# [a,b] interval of interest
# fun function
# tolerance desired accuracy
# maxIterations : maximum number of iteration
# returns:
# zero, residual, number of iterations
if (a >= b) :
print(' b must be greater than a (b > a)')
return 0,0,0
# what we consider as "zero"
eps = 1e-12
# evaluate f at the endpoints
fa = fun(a)
fb = fun(b)
if abs(fa) < eps : # a is the solution
zero = a
esterr = fa
k = 0
return zero, esterr, k
if abs(fb) < eps : # b is the solution
zero = b
esterr = fb
k = 0
return zero, esterr, k
if fa*fb > 0 :
print(' The sign of FUN at the extrema of the interval must be different')
return 0,0,0
# We want the final error to be smaller than tol,
# i.e. k > log( (b-a)/tol ) / log(2)
nmax = int(np.ceil(np.log( (b-a)/tol ) / np.log(2)))
# but nmax shall be smaller the the nmaximum iterations asked by the user
if ( maxIterations < nmax ) :
nmax = int(round(maxIterations))
print('Warning: nmax is smaller than the minimum number of iterations necessary to reach the tolerance wished');
# vector of intermadiate approximations etc
x = np.zeros(nmax)
# initial error is the length of the interval.
esterr = (b - a)
# do not need to store all the a^k and b^k, so I call them with a new variable name:
ak = a
bk = b
# the values of f at those points are fa and fk
for k in range(nmax) :
# approximate solution is midpoint of current interval
x[k] = (ak + bk) / 2
fx = fun(x[k]);
# error estimator is the half of the previous error
esterr = esterr / 2
# if we found the solution, stop the algorithm
if np.abs(fx) < eps :
# error is zero
zero = x[k]
esterr = 0;
return zero, esterr, k
if fx*fa < 0 : # alpha is in (a,x)
bk = x[k]
elif fx*fb < 0 : # alpha is in (x,b)
ak = x[k]
else :
error('Algorithm not operating correctly')
zero = x[k];
if esterr > tol :
print('Warning: bisection stopped without converging to the desired tolerance because the maximum number of iterations was reached');
return zero, esterr, k
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, inc, 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