diff --git a/NonLinearEquationsLib.py b/NonLinearEquationsLib.py index 0569538..00a9995 100644 --- a/NonLinearEquationsLib.py +++ b/NonLinearEquationsLib.py @@ -1,229 +1,229 @@ #!/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.xlabel('f'); plt.ylabel('x'); + 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','(xk,fk)'])