文章未完

相空间的绘制

我们随机选一个方程,随机选的,不是有数学手册吗,一般来说考题不可能出数学手册上的例子

import scipy.integrate as si 
import matplotlib.pyplot as plt
import numpy as np## dx/dt = x**2-y**2+x+y
## dy/dt = x*y**2 - x**2*yf  = lambda x,y:x**2-y**2+x+y
g  = lambda x,y:x*y**2 - x**2*ydt = 0.01
time_start = 0
time_end   = 1
time = np.arange(time_start,time_end+dt,dt)def system_Euler():global fglobal gglobal timex = np.zeros(time.size)y = np.zeros(time.size)x[0],y[0] = 1,2for i in range(1,len(time)):x[i] = x[i-1]+(f(x[i-1],y[i-1]))*dty[i] = y[i-1]+(g(x[i-1],y[i-1]))*dtreturn x,ydef system_odeint(X,t=0):return np.array([X[0]**2-X[1]**2+X[0]+X[1],X[0]*X[1]**2-X[0]**2*X[1]])
system_init = np.array([1,2])####################
x,y = system_Euler()fig = plt.figure(figsize=(6,6),tight_layout=True)
fig.suptitle("Stability Analysis")ax1 = plt.subplot(221)ax1.plot(time,x, 'r-', label='$x(t)$')
ax1.plot(time,y, 'b-', label='$y(t)$')
ax1.set_title("Dynamics in time")
ax1.set_xlabel("time")
ax1.grid()
ax1.legend()ax2 = plt.subplot(222)
ax2.plot(x,y,color="blue")
ax2.set_xlabel("x")
ax2.set_ylabel("y")  
ax2.set_title("Phase space")
ax2.grid()#####################33
X,infodict = si.odeint(system_odeint,system_init,time,full_output=True)
x,y = X.Tax3 = plt.subplot(223)ax3.plot(time,x, 'r-', label='$x(t)$')
ax3.plot(time,y, 'b-', label='$y(t)$')
ax3.set_title("Dynamics in time")
ax3.set_xlabel("time")
ax3.grid()
ax3.legend()ax4 = plt.subplot(224)
ax4.plot(x,y,color="blue")
ax4.set_xlabel("x")
ax4.set_ylabel("y")  
ax4.set_title("Phase space")
ax4.grid()plt.savefig("0.jpg")
plt.pause(0.01)

一阶微分方程的稳定点

import scipy.integrate as si 
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
import copy 
## dx/dt = 2*x - x**2 - x*y
## dy/dt = x*y - yf  = lambda X:X[0]*2 - X[0]**2 - X[0]*X[1]
g  = lambda X:-X[1] + X[0]*X[1]dt = 0.01
time_start = 0
time_end   = 10
time = np.arange(time_start,time_end+dt,dt)
system_init = np.array([10,2])def system_odeint(X,t=0):global fglobal greturn np.array([f(X),g(X)])fig = plt.figure(figsize=(12,6),tight_layout=True)
fig.suptitle("$x'=x(2-x-y);y'=y(-1+x)")X,infodict = si.odeint(system_odeint,system_init,time,full_output=True)
x,y = X.Tdef find_fixed_points():global fglobal gglobal Xr = sy.symbols("r")c = sy.symbols("c")R = sy.Function("R")C = sy.Function("C")R = 2*r - r**2 -r*cC = -c + r*cREqual = sy.Eq(R,0)CEqual = sy.Eq(C,0)return sy.solve((REqual,CEqual),r,c)ax1 = plt.subplot(121)
ax1.plot(time,x, 'r-', label='$x(t)$')
ax1.plot(time,y, 'b-', label='$y(t)$')
ax1.set_title("Dynamics in time")
ax1.set_xlabel("time")
ax1.grid()
ax1.legend()ax2 = plt.subplot(122)
ax2.plot(x,y,color="blue")
ax2.set_xlabel("x")
ax2.set_ylabel("y")  
ax2.set_title("Phase space")
ax2.grid()
fixed_points = find_fixed_points()
for point in fixed_points:ax2.scatter(point[0],point[1],s=20,label="fixed points")
ax2.legend()plt.savefig("0.jpg")
plt.pause(0.01)

微分方程稳定性的理论分析

• 李雅普诺夫先生，一个很棒的人

一阶自洽微分方程平衡点稳定性的结论

设有微分方程,若等号右端不显含自变量t,则称之自治方程。

代数方程的实根称为方程的平衡点（奇点，奇解）

二阶自洽微分方程平衡点稳定性的结论

一个自洽的二阶微分方程可表示为两个一阶方程组成的方程组

代数方程组的实根是方程的平衡点

一阶自洽线性常系数方程组

对于，系数矩阵

显然，方程组有唯一平衡点，假定

特征方程：

特征根：

一个典型的非线性振动系统

显然他可以转换成一个近可积哈密顿系统

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint as sio
from numpy.linalg import solve as sol
import random
random.seed(0)dt = 0.001
time_start = 0
time_end = 2
time = np.arange(time_start,time_end+dt,dt)def RK4(diff,time):global dtglobal init X = np.zeros((time.size,len(init)))X[0] = initfor i in range(time.size-1):s0 = diff(time[i],X[i])s1 = diff(time[i]+dt/2,X[i]+dt*s0/2.)s2 = diff(time[i]+dt/2,X[i]+dt*s1/2.)s3 = diff(time[i+1],X[i]+dt*s2)X[i+1] = X[i] + dt*(s0+2*(s1+s2)+s3)/6.return time,Xdef diff(time,Xi):xi,yi = Xiepision = 0.1global omegadiffxi = yidiffyi = xi**2 - xi + epision*np.cos(omega*time)return np.array([diffxi,diffyi])omega = 0.4
createvar = locals()
for i in range(20):init = np.array([1+(i>0)*random.uniform(0,0.01),1+(i>0)*random.uniform(0,0.01)])createvar["t"+str(i)],createvar["X"+str(i)] = RK4(diff,time)for i in range(20):  plt.plot(eval("X"+str(i))[:,0],eval("X"+str(i))[:,1])#,label="omega:"+str(round(omega,2)))plt.legend()plt.pause(0.01)