From 31b942368c941dd9a41f25d56e00357b9569b19c Mon Sep 17 00:00:00 2001 From: GUILHERME_MARTINS Date: Fri, 11 Jan 2019 21:16:04 -0200 Subject: [PATCH] Add bisection_method.py and Dynamic.py --- Bisection_method.py | 226 ++++++++++++++++++++++++++++++++++++++++++++ Dynamic.py | 111 ++++++++++++++++++++++ README.md | 19 +++- 3 files changed, 355 insertions(+), 1 deletion(-) create mode 100644 Bisection_method.py create mode 100644 Dynamic.py diff --git a/Bisection_method.py b/Bisection_method.py new file mode 100644 index 0000000..c0f8636 --- /dev/null +++ b/Bisection_method.py @@ -0,0 +1,226 @@ + +import numpy as np +import matplotlib.pyplot as plt + + +def bisection_2_order(Coefi,error_want = 0.00000000001,i_max=100000,mod=False): + + #Resolve second order polynomial using bisection method + #Coefi = array with coefficients. + #error_want = int. + #mod = True to return only the module of root (in case of complex number) + + + + a,b,c = Coefi + + def function(x,a=a,b=b,c=c): + f = a*x**2 + b*x + c + return f + + def sign(x,a=a,b=b,c=c): + s = function(x) + + if s < 0.0: + return (-1) + + if s > 0.0: + return (1) + + if s == 0: + return (0) + + if len(Coefi) != 3.0: + return ['ERROR','ERROR'] + + + Delta = b**2 - 4*c*a + Delta_p = round(Delta,5) + + if Delta_p < 0.0 and mod==True: + return [c,c] + if Delta_p< 0.0 and mod==False: + return ["not root","not root"] + + if Delta_p == 0: + x_1 = -b/(2*a) + x_2 = -b/(2*a) + return [x_1,x_2] + + + x_start_1 = float(-b/(2*a)) + x_start_2 = x_start_1 - 5 + + while sign(x_start_1) == sign(x_start_2): + x_start_2 = x_start_2 - 2 + + + x_start_1 = x_start_2 + 5 + media = float((x_start_2 + x_start_1)/2) + error = function(media) + i = 0 + continue_calcule = True + + + + while continue_calcule: + + + if error < error_want and error > -error_want : + x_1 = media + x_2 = c/(a*x_1) + res = False + return (x_1,x_2) + break + + if sign(x=media) != sign(x=x_start_1): + x_start_2 = media + + else: + x_start_1 = media + + media = float((x_start_2 + x_start_1)/2) + error = function(media) + + i = i + 1 + if i > i_max: + + return ["not value","not value"] + +def bisection_3_order_module(Coefi,error_want = 0.000000000001,i_max=100000,first_step=0.000001,mod=False): + + #Resolve third order polynomial using bisection method + #Coefi = array with coefficients. + #error_want = int. + #mod = True to return only the module of root (in case of complex number) + + + a,b,c,d = Coefi + a = float(a) + b = float(b) + c = float(c) + d = float(d) + def function(x,a=a,b=b,c=c,d=d): + f = a*x**3 + b*x**2 + c*x + d + return f + + def sign(x,a=a,b=b,c=c,d=d): + s = function(x) + + if s < 0.0: + return (-1) + + if s > 0.0: + return (1) + + if s == 0: + return (0) + + if len(Coefi) != 4: + return 0 + + Delta = 18*a*b*c*d - 4*d*b**3 + (b**2)*(c**2) - 4*a*c**3 - 27*(a**2)*d**2 + + + + + x_start_1,x_start_11 = bisection_2_order(Coefi=[3*a,2*b,c],mod=mod) + x_start_2_1 = x_start_1 - 1 + x_start_2_2 = x_start_1 + 1 + + + while sign(x_start_1) == sign(x_start_2_1) and sign(x_start_1) == sign(x_start_2_2): + x_start_2_1 = x_start_2_1 - first_step + x_start_2_2 = x_start_2_2 + first_step + + if sign(x_start_1) != sign(x_start_2_1): + x_start_2 = x_start_2_1 + else: + x_start_2 = x_start_2_2 + + + + media = float((x_start_2 + x_start_2 - 1)/2) + error = function(media) + i = 0 + continue_calcule = True + + + + while continue_calcule: + + + if error < error_want and error > -error_want : + x_1 = media + continue_calcule = False + break + + if sign(x=media) != sign(x=x_start_1): + x_start_2 = media + + else: + x_start_1 = media + + media = float((x_start_2 + x_start_1)/2) + error = function(media) + + i = i + 1 + if i > i_max: + print(i_max) + x_1 = media + continue_calcule=False + break + + + + if Delta < 0: + D = -d/x_1 + return [x_1,D,D] + else: + a_2 = float(a) + c_2 = float(-d/x_1) + b_2 = float(b + x_1*a) + + Coe = [a_2,b_2,c_2] + + x_2,x_3 = bisection_2_order(Coefi = Coe) + + return [x_1,x_2,x_3] + +def matrix_dete_poli(Matrix): + + #Return the coefficient of characteristic polynomial (3X3 Matrix) + + A_1,B_1,C_1 = Matrix[0] + A_2,B_2,C_2 = Matrix[1] + A_3,B_3,C_3 = Matrix[2] + a = -1 + b = A_1 + B_2 + C_3 + c = -1*C_3*(A_1 + B_2) -1*A_1*B_2 + A_3*C_1 + B_3*C_2 + A_2*B_1 + d = A_1*B_2*C_3 + B_1*C_2*A_3 + C_1*A_2*B_3 -1*B_2*A_3*C_1 -1*B_3*C_2*A_1 -1*C_3*A_2*B_1 + return [a,b,c,d] + +def eigvalues_bisection_method(Matrix,mod = False): + + #Calculate eigvalues of Matrix 3x3 . In caso complex number, need change mod to return module of complex number + + Coefi = matrix_dete_poli(Matrix) + x_1,x_2,x_3 = bisection_3_order_module(Coefi=Coefi) + return [x_1,x_2,x_3] + +def main_example_bisection_second_order(): + + #Resolvi 3x^2 + 6x - 18 + Coefi = [3,6,-18] + x_1,x_2 = bisection_2_order(Coefi=Coefi) + print(x_1,x_2) + +def main_example_bisection_third_order(): + + #Resolvi x^3 + 60x^2 - 316x - 3840 + Coefi = [1,60,-316,-3840] + x_1,x_2,x_3 = bisection_3_order_module(Coefi=Coefi) + print(x_1,x_2,x_3) + + + diff --git a/Dynamic.py b/Dynamic.py new file mode 100644 index 0000000..1b24aa6 --- /dev/null +++ b/Dynamic.py @@ -0,0 +1,111 @@ +from __future__ import division +import matplotlib.pyplot as plt +import numpy as np + + + +def Solution_discrete_system(Function,Initial_position,Parameters,Steps): + + #Return the value of system in each step and time + #Return Array NxM. N = dimensional of system. M = number of steps. + #Function = Array + #Initial_position = Array + #Parameters = Array of array. Each Parameters'function need to be array. In funciton, need separete with parameters will be used + #Steps = Int + + + + if len(Function) != len(Initial_position): + print ("No corret") + return 0 + + Position = [[]] + Future = [] + + for i in Initial_position: + Position[0].append(i) + Future.append(0) + + + + for i in range(Steps): + Position.append([]) + for j in range(len(Function)): + Future[j] = Function[j](Position[i],Parameters) + Position[i + 1].append(Future[j]) + + + + Time = range(Steps + 1) + + return Time, Position + +def plot_all_X(X,Y,labels,Title,xlim=0,ylim=0,lim=False,loc=4,save=False,missplot=False,number_plot_miss=0): + + #Plot all the variables with labels + #missplot = ignore one of variable + #number_plot_miss = the number of variable will be ignore + + + + + Separete_Y = [] + for i in range(len(Y[0])): + y = [ j[i] for j in Y ] + Separete_Y.append(y) + + + if not missplot: + for y,la in zip(Separete_Y,labels): + plt.plot(X,y,label=la) + + else: + for n_y,n_la in zip(range(len(Separete_Y)),range(len(labels))): + if n_y not in number_plot_miss: + plt.plot(X,Separete_Y[n_y],label=labels[n_la]) + try: + if not lim: + if len(xlim) == 2: + plt.xlim(left=xlim[0]) + plt.xlim(right=xlim[1]) + else: + plt.xlim(left=xlim[0]) + + if len(ylim) == 2: + plt.xlim(bottom=ylim[0]) + plt.xlim(top=ylim[1]) + else: + plt.xlim(bottom=ylim[0]) + except: + pass + + + + plt.title(Title) + plt.legend() + + if save: + plt.savefig(Title + ".png") + else: + plt.show() + + +def main_example(): + + #Define the function + def Logistic_equation(X,Parameters): + x = X[0] + R = Parameters[0] + return R*x*(1-x) + + #Create the inputs + Function = [Logistic_equation] + Initial_position = [0.1] + Parameters = [2] + Steps = 1000 + labels = [r'$x_t$'] + Title = "logistic_map" + + Time, Solution = Solution_discrete_system(Function,Initial_position,Parameters,Steps) + plot_all_X(Time,Solution,labels,Title) + diff --git a/README.md b/README.md index 69c8db7..8f7848b 100644 --- a/README.md +++ b/README.md @@ -1 +1,18 @@ -Plot cobweb, orbit diagram and Lyapunov exponents +Essencial function for working with dynamic system + +3D_plot: + cobweb() = return cobweb plot + orbit_diagram() = return orbit plot + Lyap() = return Lyapunov exponent plot + + +Bisection_method: + bisection_2_order() = return root of second order polynomial + bisection_3_order() = return root of third order polynomial + eigvalues_bisection_method() = return eigvalues of 3x3 matrix + +Dynamic: + Solution_discrete_system() = solve the dynamic + plot_all_X() = plot all variables + +