Background:
It consists of two parts, one is to show how to write code to solve ordinary differential equation, the other is to show how python calls the function to solve ordinary differential equation.
The following equations give the velocity of Lorentz attractor in three directions, and solve the motion trace
Software:
python3.5.2
Section 1: (div)
Code:
# -*- coding: utf8 -*- import numpy as np """ //Moving equation: t Location of time P(x,y,z) steps: dt Size sets: Related parameters """ def move(P, steps, sets): x, y, z = P sgima, rho, beta = sets # Speed approximation in all directions dx = sgima * (y - x) dy = x * (rho - z) - y dz = x * y - beta * z return [x + dx * steps, y + dy * steps, z + dz * steps] # Set the sets parameter sets = [10., 28., 3.] t = np.arange(0, 30, 0.01) # Location 1: P0 = [0., 1., 0.] P = P0 d = [] for v in t: P = move(P, 0.01, sets) d.append(P) dnp = np.array(d) # Location 2: P02 = [0., 1.01, 0.] P = P02 d = [] for v in t: P = move(P, 0.01, sets) d.append(P) dnp2 = np.array(d) """ //Drawing """ from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt fig = plt.figure() ax = Axes3D(fig) ax.plot(dnp[:, 0], dnp[:, 1], dnp[:, 2]) ax.plot(dnp2[:, 0], dnp2[:, 1], dnp2[:, 2]) plt.show()
Result:
Part 2:
Code:
# -*- coding: utf-8 -*- import numpy as np from scipy.integrate import odeint """ //The ordinary differential equation is defined, and the derivative of each direction, i.e. velocity, is given """ def dmove(Point,t,sets): """ p: Position vector sets: Other parameters """ p,r,b = sets x,y,z = Point return np.array([p*(y-x),x*(r-z),x*y-b*z]) t = np.arange(0,30,0.01) #Call odeint to solve dmove with two different initial values P1 = odeint(dmove,(0.,1.,0.),t,args = ([10.,28.,3.],)) #(0., 1., 0.) is the initial value of point #([10,28,3.],) parameters after point,t are given in the form of ancestor P2 = odeint(dmove,(0.,1.01,0.),t,args = ([10.,28.,3.],)) """ //Draw the curve of 3D space """ from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt fig = plt.figure() ax = Axes3D(fig) ax.plot(P1[:,0],P1[:,1],P1[:,2]) ax.plot(P2[:,0],P2[:,1],P2[:,2]) plt.show()
Result:
Copyright notice: This is the original article of the blogger. It can't be reproduced without the permission of the blogger. https://blog.csdn.net/u011702002/article/details/78118857