A simple example of solving ordinary differential equations with python and an example of solving ordinary differential equations with odeint

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

Tags: Python

Posted on Thu, 13 Feb 2020 14:05:58 -0500 by spxmgb