#! /usr/bin/env python3 # def lorenz_ode_sensitivity_test ( ): #*****************************************************************************80 # ## lorenz_ode_sensitivity_test() shows sensivitity of Lorenz ODE solutions. # # Discussion: # # This program demonstrates that small changes in the initial # condition for the Lorenz equations can result in enormous changes # in the subsequent solution trajectories. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 January 2020 # # Author: # # John Cook. # Modifications by John Burkardt. # # Reference: # # John Cook, # A different view of the Lorenz system, # https://www.johndcook.com/blog/ # 26 January 2020. # from scipy import linspace from scipy.integrate import solve_ivp import matplotlib.pyplot as plt import platform print ( '' ) print ( 'lorenz_ode_sensitivity_test:' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' Demonstrate sensitive dependence on initial data for the Lorenz system.' ) print ( ' This program is based on the work of John D Cook.' ) a = 0.0 b = 40.0 t = linspace ( a, b, 4000 ) sol1 = solve_ivp ( lorenz_deriv, [a, b], [1,1,1], t_eval = t ) sol2 = solve_ivp ( lorenz_deriv, [a, b], [1,1,1.00001], t_eval = t ) # # Plot phase portrait (x1,z1). # plt.plot ( sol1.y[0], sol1.y[2] ) plt.xlabel ( "\$x\$" ) plt.ylabel ( "\$z\$" ) filename = "lorenz_ode_sensitivity_test_x1z1.png" print ( ' Graphics saved as "', filename, '".' ) plt.savefig ( filename ) plt.show ( block = False ) plt.close ( ) # # Plot evolution (t,x1) and error evolution (t,x1-x2). # plt.subplot ( 211 ) plt.plot ( sol1.t, sol1.y[0] ) plt.grid ( True ) plt.xlabel ("\$t\$" ) plt.ylabel ("\$x_1(t)\$" ) plt.subplot ( 212 ) plt.plot ( sol1.t, sol1.y[0] - sol2.y[0] ) plt.grid ( True ) plt.xlabel ( "\$t\$" ) plt.ylabel ( "\$x_1(t) - x_2(t)\$" ) filename = "lorenz_ode_sensitivity_test_xdif.png" print ( ' Graphics saved as "', filename, '".' ) plt.savefig ( filename ) plt.show ( block = False ) plt.close ( ) # # Terminate. # print ( '' ) print ( 'lorenz_ode_sensitivity_test:' ) print ( ' Normal end of execution.' ) return def lorenz_deriv ( t, xyz ): #*****************************************************************************80 # ## lorenz_deriv() evaluates the Lorenz ODE right hand side. # # Discussion: # # The values of parameters sigma, rho and beta correspond to those used by # Lorenz in his original report. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 26 January 2020 # # Author: # # John Burkardt # # Input: # # real T: the current time. # # real XYZ[3]: the current values of X, Y, Z. # # Output: # # real DXYZDT[3]: the current values of X', Y', Z'. # x, y, z = xyz sigma = 10.0 rho = 28.0 beta = 8.0 / 3.0 dxyzdt = [ \ sigma * ( y - x ), \ x * ( rho - z ) - y, \ x * y - beta * z ] return dxyzdt def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) lorenz_ode_sensitivity_test ( ) timestamp ( )