python solution of ordinary differential equations (Systems)

Posted by bam2550 on Wed, 12 Jan 2022 20:25:19 +0100

At present, I am in the third year of junior high school. My ability is limited. If there are deficiencies, I hope to give more advice.

I saw one a week ago video So I want to use python to solve this problem.

○ analysis

Suppose there is a charged particle q in the plane with a mass of m. There is a uniform magnetic field B in space, which is perpendicular to the plane and inward, that is, a gravity field along the negative half axis of the z axis and a gravity field along the negative half axis of the y axis. Charged particles are released from point O in the magnetic field.

The equations of motion of particles can be listed directly Decompose the equation into x and y directions The solution of the equations can be obtained simultaneously.

1, dsolve method in Symphony

#Import
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
init_printing()

First declare the symbols x,y,q,m,B,g

#parameter
q,m,B,t,g = symbols('q m B t g',real=True,nonzero=True,positive=True)
#function
x,y = symbols('x y',real=True,nonzero=True,positive=True,cls=Function)

Then express the differential equation

eq1 = Eq(x(t).diff(t,t),-q*B*y(t).diff(t)/m)
eq2 = Eq(y(t).diff(t,t),-g+q*B*x(t).diff(t)/m)
sol = dsolve([eq1,eq2])

Now print out sol (better with Jupiter notebook)

sol Obviously, this formula is very complicated and simplified by trigsimp()

x = trigsimp(sol.rhs)
y = trigsimp(sol.rhs)
x,y Then, you can calculate the integral constant inside

#Define integral variables to avoid error reporting
var('C1 C2 C3 C4')
#x(0)=0
e1 = Eq(x.subs({t:0}),0)
#x'(0)=0, put subs after diff
e2 = Eq(x.diff(t).subs({t:0}),0)
#y(0)=0
e3 = Eq(y.subs({t:0}),0)
#y'(0)=0
e4 = Eq(y.diff(t).subs({t:0}),0)
l = solve([e1,e2,e3,e4],[C1,C2,C3,C4])
l Then, by substituting the integral constants into x and y, we get the final form of the solution

x = x.subs(l)
y = y.subs(l)
x,y Of course, this solution can also be written in this form Use plt to draw pictures

ts = np.linspace(0,15,1000)
consts = {q:1,B:1,g:9.8,m:1}
fx = lambdify(t,x.subs(consts),'numpy')
fy = lambdify(t,y.subs(consts),'numpy')
plt.plot(fx(ts),fy(ts))
plt.grid()
plt.show() But sympy has a disadvantage. When the differential equation is very complex, it will strike directly.

So another new method solved the problem for us

2, SciPy odeint method in integrate

#Import
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from math import e

We first create a function so that it can represent this differential equation.

For what form this function should have, start with the first-order differential equation

For example, we want to solve the following equation Its solution is obtained by simply separating variables Let's solve it with odeint, let y(1)=1

#first order
def f(x,y):
#The first derivative is expressed as a function of y and x
dydx = -y/x
return dydx
#initial condition
y0 = 1
#The initial value condition is taken at the lower limit of the independent variable x. For example, y(1) will generate a range starting from 1
x = np.linspace(1,5,100)
plt.xlim(0,5)
plt.ylim(0,5)
#For the odeint() method, the tfirst attribute refers to that the first parameter of function f is an argument
sol = odeint(f,y0,x,tfirst=True)
plt.plot(x,sol[:,0])
plt.grid()
plt.show() Similarly, for the following first-order differential equations At this point, let's set the vector u=(x,y) (column vector) and turn the system of equations into In fact, any first-order differential equation can be written in this form f means Then the differential equation can be reduced to And f(t,u) is the function we're looking for

Support vector input in odeint, so this function can be constructed like this

def f(t,u):
#x1,x2...xn = u
x,y = u
#dx1dt=...,dx2dt=...,dxndt=...
dxdt = 3*x-x*y
dydt = 2*x-y+e**t
#dudt = [dx1dt,dx2dt,...dxndt]
dudt = [dxdt,dydt]
return dudt
#x1(0)=...,x2(0)=...,xn(0)=...
u0=[0,0]
t = np.linspace(0,10,100)
sol = odeint(f,u0,t,tfirst=True)

Therefore, for second-order ordinary differential equations We first solve the equation in this form At this point, this is a system of differential equations about dy/dx and y

Let u=(y,dy/dx), we have

def f(x,u):
y,dydx = u
d2ydx2 = np.exp(x)-4*dydx-4*y
dudx = [dydx,d2ydx2]
return dudx
u0=[0,0]
x = np.linspace(0,10,100)
sol = odeint(f,u0,x,tfirst=True)

The same is true for higher order differential equations

At this time, we look back at the original problem and can easily get

def f(t,r,k,g):
#k = B*q/m
x,y,dxdt,dydt = r
d2xdt2 = -k*dydt
d2tdt2 = -g + k*dxdt
return [dxdt,dydt,d2xdt2,d2ydt2]

Define some constants

t = np.linspace(0,15,1000)
k = 1
g = 9.8

Using odeint method

r0 = [0,0,0,0]
sol = odeint(f,r0,t,tfirst = True,args=(k,g))

Draw the image

plt.plot(sol[:,0],sol[:,1])
plt.grid()
plt.show() odeint is also very accurate in solving this equation, which is almost the same as sympy.

Finally, we can let him achieve animation effects

def update_points(num):
point_ani.set_data(s[:,0][num], s[:,1][num])
return point_ani,
plt.xlabel('X')
plt.ylabel('Y')
fig,ax = plt.subplots()
plt.plot(s[:,0],s[:,1])
point_ani, = plt.plot(s[:,0], s[:,1], "ro")
plt.grid(ls="--")
# Generate animation
ani = animation.FuncAnimation(fig, update_points,range(1, len(t)), interval=5, blit=True)
plt.show() 