理论力学专题:哈密顿力学(II)----泊松括号

• 在经典力学中，泊松括号的定义是两个动态变量之间的反对称可交换二次式。泊松括号通常被用来表示物理量之间的相互作用关系，如力和速度之间的相互作用。它在哈密顿力学中有着重要的应用，可以用于描述物理系统的演化、约束和守恒规律等方面。在量子力学中，泊松括号被泛化为泊松括号算符，用于描述量子力学中的相互作用关系。
• 正则变量是理论力学中的一种量度方法，它的主要性质包括：
• 独立的：正则变量中的每个变量是相互独立
• 完备的：正则变量组成的集合是完备

泊松括号和泊松定理

• A B均是正则变量

•  泊松括号的基本性质

泊松括号改写正则方程

• 用泊松括号改写正则方程

•  用泊松括号表示任意函数的微商

泊松定理

泊松定理的内容

如果已知两个运动积分：

 存在第三个运动积分：

泊松定理的局限性

利用泊松定理，可以帮我们寻找新的运动积分。但实际上，它常常只能给出原有积分的线性组合或者恒等式，不能提供新的积分。例如：

• 理论力学中的泊松括号是一种运算符，用于描述相空间量的变化率。具体来说，泊松括号是一种双线性、反对称、满足雅可比恒等式的运算符，用于描述两个相空间量之间的变化率。泊松括号在描述诸如质点、力学系统、波动方程等物理现象中非常有用，并且也在数学上有广泛的应用。在哈密顿力学中，泊松括号经常被用于描述相空间函数的运动方程，它的重要性不言而喻。

Python 中的泊松括号

import sympy
from sympy import symbols
from sympy import diff

def poisson_bracket(f, g, x, p):
    

    dfdx = diff(f, x)
    dgdp = diff(g, p)
    dfdp = diff(f, p)
    dgdx = diff(g, x)
    return dfdx * dgdp - dfdp * dgdx



x, p = symbols('x p')
f = x ** 2 + p ** 2
g = x * p
result = poisson_bracket(f, g, x, p)
print(result)

water 一下：

泊松

• 物理学家泊松是法国数学家、物理学家、天文学家和概率论家。他为社会提出了许多重大贡献，其中包括提出了今天称为“泊松分布”和“泊松方程”的概念。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.sparse as ss
from scipy.sparse.linalg import spsolve
 
class PDE2DModel:
    #均应该传入一个一维二元数组，表示起止值
    def __init__(self,x,y):
        assert len(x)==len(y)==2,"ERROR:UNEXPECTED SV and EV!!"
        self.x = x
        self.y = y
 
    #hx表示X上的步长
    #hy表示Y上的步长
    def space_grid(self,hx,hy):
        M = int(round((self.x[1]-self.x[0])/hx,0))
        N = int(round((self.y[1]-self.y[0])/hy,0))
        assert M==N>=3,"至少网格数是合理的"
        X = np.linspace(self.x[0],self.x[1],M+1)
        Y = np.linspace(self.y[0],self.y[1],N+1)
        return M,N,X,Y
 
    def f(self,X,Y):
        return 6*X*Y**3+6*X**3*Y+np.e**X*np.sin(Y)-np.e**X*np.sin(Y)
 
    def solution(self,X,Y):
        return np.e**X*np.sin(Y)-X**3*Y**3
 
    #左边界
    def left(self,Y):
        return np.sin(Y)
    #右边界
    def right(self,Y):
        return np.e**3*np.sin(Y)-27*Y**3
    #上边界
    def up(self,X):
        return np.sin(1)*np.e**X-X**3
    #下边界
    def down(self,X):
        return 0*X
 
#解算核心
def NDM5_2D(PDE2DModel,hx,hy):
    M,N,X0,Y0 = PDE2DModel.space_grid(hx,hy)
    Y,X = np.meshgrid(Y0,X0)
##    print("X0",X0)
##    print("Y0",Y0)
##    #数值结果保存在U中 从0到N共N+1个
##    print("M",M)
##    print("N",N)
    U = np.zeros((M+1,N+1))
    U[0,:]  = PDE2DModel.left(Y0)
    U[-1,:] = PDE2DModel.right(Y0)
    U[:,0]  = PDE2DModel.down(X0)
    U[:,-1] = PDE2DModel.up(X0)
 
    D = np.diag([-1/(hy**2) for i in range(M-1)])
    C = np.zeros((N-1,N-1),dtype="float64")
    for i in range(N-1):
        C[i][i] = 2*(1/hx**2+1/hy**2)
        if i<N-2:
            C[i][i+1] = -1/hx**2
            C[i+1][i] = -1/hx**2
 
 
    u0 = np.array([[PDE2DModel.down(X0[i])] for i in range(1,M)])
    un = np.array([[PDE2DModel.up(X0[i])] for i in range(1,M)])
    
    F = np.zeros((M-1)*(N-1)) 
    for j in range(1,N):
        #for i in range(1,M):
        F[(N-1)*(j-1):(N-1)*(j)] = PDE2DModel.f(X0[1:M],np.array([Y0[j] for i in range(N-1)]))
 
        F[(N-1)*(j-1)] += PDE2DModel.left(Y0[j])/hx**2
        F[(N-1)*(j)-1] += PDE2DModel.right(Y0[j])/hx**2
 
    F[:N-1] -= np.dot(D,u0).T[0]
    F[(M-1)*(N-2):] -= np.dot(D,un).T[0]
    F = np.mat(F).T
##    print(F)
 
    Dnew = np.zeros(((M-1)*(N-1),(N-1)*(M-1)))
    for i in range((N-1)*(N-1)):
        Dnew[i][i] = 2*(1/hx**2+1/hy**2)
 
        if i<(N-2)*(N-1):
            Dnew[i][i+N-1] = -1/hy**2
            Dnew[i+N-1][i] = -1/hy**2
 
    for i in range(N-1):
        for j in range(N-2):
            Dnew[(N-1)*i+j][(N-1)*i+j+1] = -1/hx**2
            Dnew[(N-1)*i+j+1][(N-1)*i+j] = -1/hx**2
 
    print("差分方程构造完成！解算开始！")  
    Unew = np.linalg.solve(Dnew,F)
    #print(Unew)
    U[1:-1,1:-1] = Unew[:,0].reshape((N-1,N-1)).T
    return X,Y,U
 
#数据可视化
def Visualized():
    x = np.array([0,3])
    y = np.array([0,1])
 
    pde = PDE2DModel(x,y)
    X,Y,U = NDM5_2D(pde,0.03,0.01)
    u = pde.solution(X,Y)
 
    print("解算完成！绘图已开始！")
    plt.figure(figsize=(15,5))
    ax1 = plt.subplot(131,projection="3d")
    ax2 = plt.subplot(132,projection="3d")
    ax3 = plt.subplot(133,projection="3d")
 
    ax1.set_title("Numeric Solution")
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.plot_surface(X,Y,U,cmap="gist_ncar")
 
    ax2.set_title("Exact Solution")
    ax2.set_xlabel("x")
    ax2.set_ylabel("y")
    ax2.plot_surface(X,Y,u,cmap="gist_ncar")
 
    e = np.abs(U-u)
    ax3.set_title("Error")
    ax3.set_xlabel("x")
    ax3.set_ylabel("y")
    ax3.plot_surface(X,Y,e,cmap="gist_ncar")
 
    plt.show()
    return U,u,X,Y
 
U,u,X,Y = Visualized()