精华内容
下载资源
问答
  • 首先介绍一个特殊的方阵J∈R2n×2nJ \in R^{2n \times 2n}J∈R2n×2n,其定义如下: J=[OnIn−InOn]J=\begin{bmatrix} O_n & I_n\\ -I_n & O_n \end{bmatrix}J=[On​−In​​In​On​​] 其中On∈Rn×nO_n\...

    首先介绍一个特殊的方阵 J ∈ R 2 n × 2 n J \in R^{2n \times 2n} JR2n×2n,其定义如下:
    J = [ O n I n − I n O n ] J=\begin{bmatrix} O_n & I_n\\ -I_n & O_n \end{bmatrix} J=[OnInInOn]
    其中 O n ∈ R n × n O_n\in R^{n\times n} OnRn×n是零矩阵; I n ∈ R n × n I_n \in R^{n \times n} InRn×n是单位矩阵。

    那么这个矩阵具有下列的性质:
    J T = − J J^T = -J JT=J
    J − 1 = J T J^{-1} = J^T J1=JT
    J T J = I 2 n J^TJ = I_{2n} JTJ=I2n
    J T J T = − I 2 n J^TJ^T = -I_{2n} JTJT=I2n
    J 2 = − I 2 n J^2 = -I_{2n} J2=I2n
    d e t J = ± 1 det J = \pm 1 detJ=±1

    哈密顿矩阵(Hamiltonian matrix)定义

    一个矩阵 A ∈ R 2 n × 2 n A\in R^{2n \times 2n} AR2n×2n叫哈密顿矩阵如果 J A JA JA 是对称的,也就是说:
    J A = ( J A ) T ⇒ A T J + J A = 0 JA = (JA)^T \Rightarrow A^TJ+JA = 0 JA=(JA)TATJ+JA=0,其中 J J J就是上面介绍的特殊矩阵。
    H = { A ∈ R 2 n × 2 n ∣ A T J + J A = 0 } \mathcal{H} =\{A \in R^{2n\times 2n}|A^TJ+JA = 0\} H={AR2n×2nATJ+JA=0} 2 n × 2 n 2n\times 2n 2n×2n的哈密顿矩阵的集合。

    以下三条性质是等价的:

    1. A A A是一个哈密顿矩阵
    2. A = J S A=JS A=JS,其中 S = S T S=S^T S=ST
    3. (JA)^T = JA

    定理:
    A , B ∈ H n A,B \in \mathcal{H}^n A,BHn,那么下面定理是对的:

    1. A + B ∈ H n A+B \in \mathcal{H} ^n A+BHn
    2. α A ∈ H n \alpha A \in \mathcal{H}^n αAHn
    3. [ A , B ] ∈ H n [A,B] \in \mathcal{H}^n [A,B]Hn,其中 [ A , B ] = d e f A B − B A [A,B]\overset{\underset{\mathrm{def}}{}}{=}AB-BA [A,B]=defABBA

    结果: ( H , [ ⋅ , ⋅ ] ) (\mathcal{H}, [\cdot, \cdot]) (H,[,])是个李代数(Lie algebra)。

    定理

    A ∈ H n A\in \mathcal{H}^n AHn p A ( x ) p_A(x) pA(x)为矩阵 A A A的特征多项式,那么:

    • p A ( x ) = p A ( − x ) p_A(x) = p_A(-x) pA(x)=pA(x)
    • 如果 p A ( c ) = 0 p_A(c) = 0 pA(c)=0,那么 p A ( − c ) = p A ( c ˉ ) = p A ( − c ˉ ) = 0 p_A(-c) = p_A(\bar c)=p_A(-\bar c)=0 pA(c)=pA(cˉ)=pA(cˉ)=0,其中 c ∈ R c\in R cR

    代数Riccati Equation

    首先定义不变子空间(invariant subspace)
    由向量 v 1 , , . . . v k v_1,,...v_k v1,,...vk张成的线性空间 ν \nu ν叫做矩阵 A A A的不变子空间如果对于任意 v ∈ ν , A v ∈ ν v \in \nu, Av \in \nu vν,Avν
    考虑下面的代数Riccati Equation (ARE):
    0 = R ( x ) = F + A T X + X A + X A − X G X ( 1 ) 0=R(x)=F+A^TX+XA+XA-XGX (1) 0=R(x)=F+ATX+XA+XAXGX1
    其中 A , F , G , X ∈ R n × n A,F,G,X\in R^{n \times n} A,F,G,XRn×n F , G F,G F,G是对称矩阵, X X X是未知的对称矩阵。
    定义 H = [ A G F − A T ] H =\begin{bmatrix}A & G\\ F& -A^T \end{bmatrix} H=[AFGAT],令 [ U V ] \begin{bmatrix} U\\ V \end{bmatrix} [UV]
    H H H的不变子空间的一个向量,
    也就是说
    [ A G F − A T ] [ U V ] = [ U V ] Z ,        ( 2 ) \begin{bmatrix} A & G\\ F & -A^T \end{bmatrix}\begin{bmatrix} U\\ V \end{bmatrix}=\begin{bmatrix} U\\ V \end{bmatrix}Z , \; \; \; (2) [AFGAT][UV]=[UV]Z,(2)
    其中 Z ∈ R n × n , λ ( Z ) ⊂ λ ( H ) Z\in R^{n\times n}, \lambda(Z) \subset \lambda(H) ZRn×n,λ(Z)λ(H)
    假定 U U U是非奇异的,那么从 ( 2 ) (2) (2)的第一行得到:
    A U + G V = U Z ↔ U − 1 A U + U − 1 G V = Z AU+GV = UZ \leftrightarrow U^{-1}AU+U^{-1}GV=Z AU+GV=UZU1AU+U1GV=Z
    将其插入第二行得到的方程:
    F U − A T V = V Z = V U − 1 A U + V U − 1 G V FU-A^TV=VZ=VU^{-1}AU+VU^{-1}GV FUATV=VZ=VU1AU+VU1GV
    以上方程等于
    0 = F − A T V U − 1 − V U − 1 A − V U − 1 G V U − 1 0=F-A^TVU^{-1}-VU^{-1}A-VU^{-1}GVU^{-1} 0=FATVU1VU1AVU1GVU1
    X : = − V U − 1 X:= -VU^{-1} X:=VU1,我们发现 X X X正好是 ( 1 ) (1) (1)的解。

    定理

    考虑哈密顿矩阵 H H H,假设该矩阵没有特征值在虚轴上,其不变子空间是 χ − = I m [ X 1 X 2 ] ⊂ R 2 n × n \chi_{-} = Im \begin{bmatrix} X_1\\ X_2 \end{bmatrix} \subset R^{2n \times n} χ=Im[X1X2]R2n×n
    这是一个由稳定特征值对应的特征向量张成的空间。
    如果 X 1 X_1 X1是可逆的,那么 X = − X 2 X 1 − 1 X=-X_2X_1^{-1} X=X2X11是代数Ricatti方程的解,并且 A − G X A-GX AGX是稳定的。

    展开全文
  • 量子化学哈密顿矩阵元的计算,基于上一篇文章构建的占据数矢量。

    数学推导

    本文是上一篇文章 http://blog.csdn.net/luozhen07/article/details/78701311 的继续。

    量子化学计算的核心过程是构建体系的哈密顿矩阵并对角化,从而得到能量(特征值)和对应的波函数(特征矢)。波函数由一个或多个占据数矢量线性组合而成:前者称为“单行列式波函数”,譬如最基础的限制性 Hartree-Fock 自洽场(Hartree-Fock SCF)方法以及密度泛函(Density Functional Theory)方法得到的波函数;后者称为“多行列式波函数”,例如多组态自洽场(Multi-Configuration SCF, MCSCF)方法以及这里要讨论的组态相互作用(Configuration Interaction, CI)方法得到的波函数。用符号表示,为

    |Ψ=kckϕk

    哈密顿矩阵元可以表示为

    HMN=ΨM|H^|ΨN=i,jcM,icN,jϕi|H^|ϕj

    即,最终需要计算两个行列式之间的哈密顿量。考虑哈密顿算符的形式
    H^=pqhpqEpq+12pqrsgpqrsepqrs

    其中,单激发算符为
    Epq=apαaqα+apβaqβ

    双激发算符为
    epqrs=apαarαasαaqα+apαarβasβaqα+apβarαasαaqβ+apβarβasβaqβ

    并且 hpq gpqrs 分别为单、双电子积分。考虑
    ϕi|ϕj=δij

    可见,只有当两个行列式的差异不超过两次激发(即,通过不超过两次电子移动,可以让左右两个行列式完全相同)时,哈密顿量才有非零值。这部分的具体讨论,可以参考下面的代码或者搜索“ Slater-Condon 规则”。

    程序实现

    这里的程序主要是用于计算行列式(或者说“占据数矢量”)之间的哈密顿矩阵元 ϕi|H^|ϕj 。在上面的数学推导中,我们知道此计算还需要单、双电子积分。在这部分代码中,我们使用了名为 Integral 的类来保存积分,并且 Integral::h1e(p,q) 函数返回 p 轨道与 q 轨道之间的单电子积分,Integral::h2e(p,q,r,s) 函数返回 pqrs 轨道之间的双电子积分。除了计算哈密顿矩阵元,其他计算主要有计算电子激发个数和激发的位置。

    /* --- file: ci.h
       --- author: LUO Zhen,
                   Institute of Theoretical and Computational Chemistry,
                   Nanjing University
    */
    
    #ifndef CONFIGURATION_INTERACTION_FUNCTIONS_HEADER
    #define CONFIGURATION_INTERACTION_FUNCTIONS_HEADER
    
    #include "slaterdet.h"
    #include "integral.h"
    #include "basicops.h"
    
    // excitations are always from the second parameter (det_j) to the first one (det_i).
    
    namespace CI
    {
        // number of excitations from det_j to det_i, alpha
        int n_excitations_a(const SlaterDet& det_i, const SlaterDet& det_j);
        // number of excitations from det_j to det_i, beta
        int n_excitations_b(const SlaterDet& det_i, const SlaterDet& det_j);
        // number of excitations from det_j to det_i, alpha + beta
        int n_excitations(const SlaterDet& det_i, const SlaterDet& det_j);
        // compute orbital indices for a single excitation: alpha, return int[2]
        int* get_single_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j);
        // compute orbital indices for a single excitation: beta, return int[2]
        int* get_single_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j);
        // compute orbital indices for a double excitation: alpha, return int[4]
        int* get_double_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j);
        // compute orbital indices for a double excitation: beta, return int[4]
        int* get_double_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j);
        // compute Hamiltonian matrix element H_ij = <det_i|H|det_j>
        double H_ij(const SlaterDet& det_i, const SlaterDet& det_j, const Integral& integrals, int norb);
    }
    
    #endif // !CONFIGURATION_INTERACTION_FUNCTIONS_HEADER

    这些函数的实现如下:

    /* --- file: ci.cpp
       --- author: LUO Zhen,
                   Institute of Theoretical and Computational Chemistry,
                   Nanjing University
    */
    
    #include "ci.h"
    
    // number of excitations from det_j to det_i, alpha
    int CI::n_excitations_a(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::n_excitations(det_i.alpha(), det_j.alpha(), SlaterDet::nset());
    }
    
    // number of excitations from det_j to det_i, beta
    int CI::n_excitations_b(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::n_excitations(det_i.beta(), det_j.beta(), SlaterDet::nset());
    }
    
    // number of excitations from det_j to det_i, alpha + beta
    int CI::n_excitations(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return n_excitations_a(det_i, det_j) + n_excitations_b(det_i, det_j);
    }
    
    // compute orbital indices for a single excitation: alpha, return int[2]
    int* CI::get_single_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::get_single_excitation(det_i.alpha(), det_j.alpha(), SlaterDet::nset());
    }
    
    // compute orbital indices for a single excitation: beta, return int[2]
    int* CI::get_single_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::get_single_excitation(det_i.beta(), det_j.beta(), SlaterDet::nset());
    }
    
    // compute orbital indices for a double excitation: alpha, return int[4]
    int* CI::get_double_excitation_a(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::get_double_excitation(det_i.alpha(), det_j.alpha(), SlaterDet::nset());
    }
    
    // compute orbital indices for a double excitation: beta, return int[4]
    int* CI::get_double_excitation_b(const SlaterDet& det_i, const SlaterDet& det_j)
    {
        return BasicOps::get_double_excitation(det_i.beta(), det_j.beta(), SlaterDet::nset());
    }
    
    // compute Hamiltonian matrix element H_ij = <det_i|H|det_j>
    double CI::H_ij(const SlaterDet& det_i, const SlaterDet& det_j, const Integral& integrals, int norb)
    {
        double h1 = 0.0; // single particle
        double h2 = 0.0; // double particle
        int n_excit_a = CI::n_excitations_a(det_i, det_j);
        int n_excit_b = CI::n_excitations_b(det_i, det_j);
        if(n_excit_a + n_excit_b < 3) // ignore >= 3 excitations
        {
            int nelec_a = det_i.nelec_a();
            int nelec_b = det_i.nelec_b();
            if(n_excit_a == 0 && n_excit_b == 0) // i == j, diagonal term
            {
                int* occ_list_a = det_i.occ_list_a(norb);
                int* occ_list_b = det_i.occ_list_b(norb);
                // single particle
                // h1 = SUM_(p in occ_a){h1e(p,p)} + SUM_(p in occ_b){h1e(p,p)}
                for(int idx1 = 0; idx1 < nelec_a; ++idx1) // alpha
                {
                    int p = occ_list_a[idx1];
                    h1 += integrals.h1e(p, p);
                }
                for(int idx1 = 0; idx1 < nelec_b; ++idx1) // beta
                {
                    int p = occ_list_b[idx1];
                    h1 += integrals.h1e(p, p);
                }
                // double particle
                // h2 = 0.5 * (SUM(p,r in occ_a){h2e(p,p,r,r) - h2e(p,r,r,p)}
                //            +SUM(p,r in occ_b){h2e(p,p,r,r) - h2e(p,r,r,p)}
                //            +SUM(p in occ_a, r in occ_b){h2e(p,p,r,r)}
                //            +SUM(p in occ_b, r in occ_a){h2e(p,p,r,r)})
                for(int idx1 = 0; idx1 < nelec_a; ++idx1)
                {
                    int p = occ_list_a[idx1];
                    // SUM(p,r in occ_a){h2e(p,p,r,r) - h2e(p,r,r,p)}
                    for(int idx2 = 0; idx2 < nelec_a; ++idx2)
                    {
                        int r = occ_list_a[idx2];
                        h2 += (integrals.h2e(p, p, r, r) - integrals.h2e(p, r, r, p));
                    }
                    // SUM(p in occ_a, r in occ_b){h2e(p,p,r,r)}
                    for(int idx2 = 0; idx2 < nelec_b; ++idx2)
                    {
                        int r = occ_list_b[idx2];
                        h2 += integrals.h2e(p, p, r, r);
                    }
                }
                for(int idx1 = 0; idx1 < nelec_b; ++idx1)
                {
                    int p = occ_list_b[idx1];
                    // SUM(p,r in occ_b){h2e(p,p,r,r) - h2e(p,r,r,p)}
                    for(int idx2 = 0; idx2 < nelec_b; ++idx2)
                    {
                        int r = occ_list_b[idx2];
                        h2 += (integrals.h2e(p, p, r, r) - integrals.h2e(p, r, r, p));
                    }
                    // SUM(p in occ_b, r in occ_a){h2e(p,p,r,r)}
                    for(int idx2 = 0; idx2 < nelec_a; ++idx2)
                    {
                        int r = occ_list_a[idx2];
                        h2 += integrals.h2e(p, p, r, r);
                    }
                }
                h2 *= 0.5;
                delete[] occ_list_b;
                delete[] occ_list_a;
            }
            else if(n_excit_a == 1 && n_excit_b == 0) // a -> a
            {
                int* ia = CI::get_single_excitation_a(det_i, det_j);
                int* occ_list_a = det_i.occ_list_a(norb);
                int* occ_list_b = det_i.occ_list_b(norb);
                // single particle
                // h1 = sign(p,q) * h1(p,q)
                int p = ia[1];
                int q = ia[0];
                double sign = BasicOps::compute_cre_des_sign(p, q, det_j.alpha(), SlaterDet::nset());
                h1 = sign * integrals.h1e(p, q);
                // double particle
                // h2 = sign(p,q) * (SUM_(r in occ_a){h2e(p,q,r,r) - h2e(p,r,r,q)}
                //                  +SUM_(r in occ_b){h2e(p,q,r,r)})
                for(int idx1 = 0; idx1 < nelec_a; ++idx1)
                {
                    int r = occ_list_a[idx1];
                    h2 += (integrals.h2e(p, q, r, r) - integrals.h2e(p, r, r, q));
                }
                for(int idx2 = 0; idx2 < nelec_b; ++idx2)
                {
                    int r = occ_list_b[idx2];
                    h2 += integrals.h2e(p, q, r, r);
                }
                h2 *= sign;
                delete[] occ_list_b;
                delete[] occ_list_a;
                delete[] ia;
            }
            else if(n_excit_a == 0 && n_excit_b == 1) // b -> b
            {
                int* ia = CI::get_single_excitation_b(det_i, det_j);
                int* occ_list_a = det_i.occ_list_a(norb);
                int* occ_list_b = det_i.occ_list_b(norb);
                // single particle
                // h1 = sign(p,q) * h1(p,q)
                int p = ia[0];
                int q = ia[1];
                double sign = BasicOps::compute_cre_des_sign(p, q, det_j.beta(), SlaterDet::nset());
                h1 = sign * integrals.h1e(p, q);
                // double particle
                // h2 = sign(p,q) * (SUM_(r in occ_b){h2e(p,q,r,r) - h2e(p,r,r,q)}
                //                  +SUM_(r in occ_a){h2e(p,q,r,r)})
                for(int idx1 = 0; idx1 < nelec_b; ++idx1)
                {
                    int r = occ_list_b[idx1];
                    h2 += (integrals.h2e(p, q, r, r) - integrals.h2e(p, r, r, q));
                }
                for(int idx2 = 0; idx2 < nelec_a; ++idx2)
                {
                    int r = occ_list_a[idx2];
                    h2 += integrals.h2e(p, q, r, r);
                }
                h2 *= sign;
                delete[] occ_list_b;
                delete[] occ_list_a;
                delete[] ia;
            }
            else if(n_excit_a == 2 && n_excit_b == 0) // a,a -> a,a
            {
                int* ijab = CI::get_double_excitation_a(det_i, det_j);
                int p = ijab[0];
                int r = ijab[1];
                int s = ijab[2];
                int q = ijab[3];
                // make sure s < q and p < r
                if(q < s) std::swap(q, s);
                if(r < p) std::swap(p, r);
                double sign = BasicOps::compute_cre_des_sign(p, r, det_i.alpha(), SlaterDet::nset());
                sign *= BasicOps::compute_cre_des_sign(s, q, det_j.alpha(), SlaterDet::nset());
                h2 += sign * (integrals.h2e(p, s, r, q) - integrals.h2e(p, q, r, s));
                delete[] ijab;
            }
            else if(n_excit_a == 0 && n_excit_b == 2) // b,b -> b,b
            {
                int* ijab = CI::get_double_excitation_b(det_i, det_j);
                int p = ijab[0];
                int r = ijab[1];
                int s = ijab[2];
                int q = ijab[3];
                // make sure s < q and p < r
                if(q < s) std::swap(q, s);
                if(r < p) std::swap(p, r);
                double sign = BasicOps::compute_cre_des_sign(p, r, det_i.beta(), SlaterDet::nset());
                sign *= BasicOps::compute_cre_des_sign(s, q, det_j.beta(), SlaterDet::nset());
                h2 += sign * (integrals.h2e(p, s, r, q) - integrals.h2e(p, q, r, s));
                delete[] ijab;
            }
            else // a,b -> a,b
            {
                int* ia = CI::get_single_excitation_a(det_i, det_j);
                int* jb = CI::get_single_excitation_b(det_i, det_j);
                int p = ia[0];
                int q = ia[1];
                int r = jb[0];
                int s = jb[1];
                double sign = BasicOps::compute_cre_des_sign(q, p, det_j.alpha(), SlaterDet::nset());
                sign *= BasicOps::compute_cre_des_sign(s, r, det_j.beta(), SlaterDet::nset());
                h2 += (sign * integrals.h2e(p, q, r, s));
                delete[] jb;
                delete[] ia;
            }
        }
        return h1 + h2;
    }

    多组态波函数哈密顿矩阵元的计算

    由之前的数学推导可知,对于多组态波函数哈密顿矩阵元,只需要对左矢和右矢分别遍历,计算每一组行列式的矩阵元与对应系数的乘积并求加和。定义如下的行列式集合类

    class DetSet
    {
    public:
        // 省略的代码 ...
    private:
        std::vector<double> _coeffs;
        std::vector<SlaterDet> _dets;
    };

    则哈密顿矩阵元可以通过如下方式计算

    double H_ij(const DetSet& set_i, const DetSet& set_j, const Integral& integrals, int norb)
    {
        double h = 0.0;
    #pragma omp parallel for reduction(+:h)
        for(int idx1 = 0; idx1 < set_i.size(); ++idx1)
        {
            for(int idx2 = 0; idx2 < set_j.size(); ++idx2)
            {
                double coeff_product = set_i.coeff(idx1) * set_j.coeff(idx2);
                h += (coeff_product * CI::H_ij(set_i.det(idx1), set_j.det(idx2), integrals, norb));
            }
        }
        return h;
    }
    展开全文
  • 哈密顿路径

    千次阅读 2019-01-10 12:33:50
    这个“哈密瓜路径”网上查了好久没搞明白,我这个的代码无奈定义了4个节点和5条边,将就混过课程设计课。 放代码和结果从我做起!欢迎大家留言! 1. 问题描述 在图G中找出一条包含所有顶点的简单路径,该路径称为...

    这个“哈密瓜路径”网上查了好久没搞明白,我这个的代码无奈定义了4个节点和5条边,将就混过课程设计课。

    放代码和结果从我做起!欢迎大家留言!

    1.  问题描述

    在图G中找出一条包含所有顶点的简单路径,该路径称为哈密顿路径。

    2. 基本要求

    (1)图G是非完全有向图,且图G不一定存在哈密顿路径;

    (2)设计算法判断图G是否存在哈密顿路径,如果存在,输出一条哈密顿路径即可;

    (3)分析算法的时间复杂度。

    3. 问题分析

    哈密顿路径也称作哈密顿链,指在一个图中沿边访问每个顶点恰好一次的路径。寻找这样的一个路径是一个典型的NP-完全(NP-complete)问题。图中有的边可以不经过,但是不会有边被经过两次。与欧拉图的区别:欧拉图讨论的实际上是图上关于边的可行便利问题,而哈密顿图的要求与点有关。

    判定:

    一.Dirac定理(充分条件)

    设一个无向图中有N个顶点,若所有顶点的度数大于等于N/2,则哈密顿回路一定存在(N/2指的是向上取整)。

    二.基本的必要条件

    设图G=<V, E>是哈密顿图,则对于v的任意一个非空子集S,若以|S|表示S中元素的数目,G-S表示G中删除了S中的点以及这些点所关联的边后得到的子图,则W(G-S)<=|S|成立。其中W(G-S)是G-S中联通分支数。

    三.竞赛图(哈密顿通路)

    N(N>=2)阶竞赛图一点存在哈密顿通路。

    4. 概要设计

    在Dirac定理的前提下构造哈密顿回路的过程:

    1)任意找两个相邻的节点S和T,在其基础上扩展出一条尽量长的没有重复结点的路径。即如果S与结点v相邻,而且v不在路径S -> T上,则可以把该路径变成v -> S -> T,然后v成为新的S。从S和T分别向两头扩展,直到无法继续扩展为止,即所有与S或T相邻的节点都在路径S -> T上。

    2)若S与T相邻,则路径S -> T形成了一个回路。

    3) 若S与T不相邻,可以构造出来一个回路.设路径S -> T上有k+2个节点,依次为S, v1, v2, ..., vk, T.可以证明存在节点vi(i属于[1, k]),满足vi与T相邻,且vi+1与S相邻.找到这个节点vi,把原路径变成S -> vi -> T -> vi+1 ,即形成了一个回路。

    4)到此为止,已经构造出来了一个没有重复节点的的回路,如果其长度为N,则哈密顿回路就找到了。如果回路的长度小于N,由于整个图是连通的,所以在该回路上,一定存在一点与回路之外的点相邻。那么从该点处把回路断开,就变回了一条路径,同时还可以将与之相邻的点加入路径。再按照步骤1的方法尽量扩展路径,则一定有新的节点被加进来。接着回到路径2。

    证明:

      根据鸽巢定理,既然与S和T相邻的点都在路径上,它们分布的范围只有v1,v2,---,vk这k个点,k<=N-2,跟据哈密顿回路的第一个判断条件,d(S)+d(T)>=N,那么v1,v2,---,vk这k个点中一定有至少2个点同时与S和T相连,根据鸽巢定理,肯定存在一个与S相邻的点vi和一个与T相邻的点vj,满足j=i+1。

    哈密顿路径的伪代码如下:

    算法:哈密顿路径

    输入:哈密顿回路的起始点s,哈密顿回路中终点s之前的点t

    输出:最终的哈密顿回路ans[ ]

               1.初始化,令s = 1,t为s的任意一个邻接点;

                2.如果ans[]中元素的个数小于n,则从t开始向外扩展,如果有可扩展点v,放入ans[]的尾部,并且t=v,并继续扩展,如无法扩展进入步骤3;

                3.将当前得到的ans[]倒置,s和t互换,从t开始向外扩展,如果有可扩展点v,放入ans[]尾部,并且t=v,并继续扩展。如无法扩展进入步骤4;

                        3.1如果当前s和t相邻,进入步骤5。否则,遍历ans[],寻找点ans[i],使得ans[i]与t相连并且ans[i +1]与s相连,将从ans[i + 1]到t部分的ans[]倒置,t=ans[i +1],进如步骤5;

                        3.2如果当前ans[]中元素的个数等于n,算法结束,ans[]中保存了哈密顿回路(可看情况是否加入点s).否则,如果s与t连通,但是ans[]中的元素的个数小于n,则遍历ans[],寻找点ans[i],使得ans[i]与ans[]外的一点(j)相连,则令s=ans[i - 1],t = j,将ans[]中s到ans[i - 1]部分的ans[]倒置,将ans[]中的ans[i]到t的部分倒置,将点j加入到ans[]的尾部,转步骤2;

     

    时间复杂度:

      如果说每次到步骤5算一轮的话,那么由于每一轮当中至少有一个节点被加入到路径S -> T中,所以总的轮数肯定不超过n轮,所以时间复杂度为O(n^2).空间上由于边数非常多,所以采用邻接矩阵来存储比较适合。

    5. 编程结果

    #include<stdio.h>
    #include<stdlib.h>
    #define Maxsize 10
    #define MAX_NO_LIMIT 10000000
    typedef char Datatype;
    int visited[Maxsize] = { 0 };
    int s[Maxsize] = { 0 };
    int count = 0;
    typedef struct
    {
    	Datatype vertex[Maxsize];
    	int edge[Maxsize][Maxsize];
    	int vertexnum, edgenum;
    
    }graph;
    void creatgraph(graph *G, Datatype a[], int n, int e)
    {
    	int i, j, k;
    	G->vertexnum = n;
    	G->edgenum = e;
    	for (i = 0; i<G->vertexnum; i++)
    	{
    		G->vertex[i] = a[i];
    	}
    	for (i = 0; i<G->vertexnum; i++)
    	{
    		for (j = 0; j<G->vertexnum; j++)
    		{
    			G->edge[i][j] = 0;
    		}
    	}
    	for (k = 0; k<G->edgenum; k++)
    	{
    		printf("输入两个定点");
    		scanf_s("%d%d", &i, &j);
    		G->edge[i][j] = 1;
    		G->edge[j][i] = MAX_NO_LIMIT;
    	}
    
    }
    void DFS(graph *G, int v)          //count是全局变量并已初始化为0
    {
    	int j;
    	visited[v] = 1; //对访问点进行标记
    	s[count++] = v;//S表示路径,这里为开头为v
    	for (j = 0; j <G->vertexnum; j++) {
    		if (G->edge[v][j] == 1 && visited[j] == 0)
    			DFS(G, j); //这里对点进行递归
    	}
    	if (j == G->vertexnum) {   //取消回溯
    		visited[v] = 0;
    		count--;
    	}
    }
    int main()
    {
    	int i;
    	char ch[] = { 'a','b','c','d' };
    	graph MG;
    	creatgraph(&MG, ch, 4, 5);
    	DFS(&MG, 0);
    	for (i = 0; i<4; i++)
    	{
    		printf("%d",s[i]);
    	}
    	system("pause");
    	return 0;
    }

     

    展开全文
  • 哈密顿回路

    千次阅读 2019-08-13 08:23:57
    哈密顿图(哈密尔顿图)(英语:Hamiltonian graph,或Traceable graph)是一个无向图,由天文学家哈密顿提出,由指定的起点前往指定的终点,途中经过所有其他节点且只经过一次。在图论中是指含有哈密顿回路的图,...

    哈密顿图(哈密尔顿图)(英语:Hamiltonian graph,或Traceable graph)是一个无向图,由天文学家哈密顿提出,由指定的起点前往指定的终点,途中经过所有其他节点且只经过一次。在图论中是指含有哈密顿回路的图,闭合的哈密顿路径称作哈密顿回路(Hamiltonian cycle),含有图中所有顶点的路径称作哈密顿路径(Hamiltonian path)。

    由来

     

    哈密顿回路

    天文学家哈密顿(William Rowan Hamilton) 提出,在一个有多个城市的地图网络中,寻找一条从给定的起点到给定的终点沿 途恰好经过所有其他城市一次的路径。

    这个问题和著名的七桥问题的不同之处在于,过桥只需要确定起点,而不用确定终点。哈密顿问题寻找一条从给定的起点到给定的终点沿 途恰好经过所有其他城市一次的路径

    算法

    哈密顿路径问题在上世纪七十年代初,终于被证明是“NP完全”的。据说具有这样性质的问题,难于找到一个有效的算法。实际上对于某些顶点数不到100的网络,利用现有最好的算法和计算机也需要比较荒唐的时间(比如几百年)才能确定其是否存在一条这样的路径。

    从图中的任意一点出发,路途中经过图中每一个结点当且仅当一次,则成为哈密顿回路。

    要满足两个条件:

    ⒈封闭的环

    ⒉是一个连通图,且图中任意两点可达

    经过图(有向图或无向图)中所有顶点一次且仅一次的通路称为哈密顿通路。

    经过图中所有顶点一次且仅一次的回路称为哈密顿回路。

    具有哈密顿回路的图称为哈密顿图,具有哈密顿通路但不具有哈密顿回路的图称为半哈密顿图。

    平凡图是哈密顿图。

    生成下一个结点,代码如下:

    void NextValue(int k) {
    //X(l),…,X(k-1)是一条有k-l个不同结点的路径。若X(k)=0,则表示再无结
    //点可分配给X(k)。若还有与X(1),…,X(k-1)不同且与X(k-1)有边相连结的
    //结点则将其中标数最高的结点置于X(k)。若k=n,则还需与X(1)相连结。
    int n,X[n];bool Graph[n][n]; // n、X[n]定义成全局变量
    int k,j;
    while(1) {
    X[k]=(x[k] + 1) mod (m+1) //下一个结点
    if(X[k]==0) { return;}
    if(Graph[X[k-1]],X[k]) { //有边相连吗
    for(j=1;j<=k-1;++j) { //检查与前k-1个结点是否相同
    if(X[j]==X[k]) break; //有相同结点,出此循环
    };//for
    if(j==k) { //若为真,则是一个不同结点
    if((k<n) or (k==n and Graph[X[n]],[1] ) return;
    };//if
    };//if
    }//while
    }//NextValue

    使用过程NextValue和将递归回溯算法细化得到算法Hamiltonian。此算法可找出所有的哈密顿环。

     

    找所有的哈密顿环,代码如下:

    void Hamiltonian(int k) {
    //这是找出图G中所有哈密顿环的递归算法。图G用它的布尔邻接矩阵
    //Graph(1..n,1..n)表示。每个环都从结点 1开始。
    int X[n]; //X[n]定义成全局变量
    int k,n;
    while(1) { //生成X(k)的值
    NextValue(k); //下一个合法结点分配给X(k)
    if(x[k]==0) return;
    if(k==) { print(X,’1’);} //打印一个环
    else { Hamiltonian(k+1)}
    //if
    }//repeat
    }//Hamiltonian

     这个过程首先初始化邻接矩阵Graph(l..n,l..n),然后置X(2:n)=0,X(l)=l,再执行Hamiltonian(2)

     

     

    展开全文
  • Kerov哈密顿定义为一组具有Kerov函数作为常见本征函数的通勤算子。 在Macdonald多项式的特殊情况下,...这些哈密顿量的问题是它们 借助Kostka矩阵(而不是对其进行定义)来构造,因此功能不如Ruijsenaars强大。
  • 哈密顿问题

    2018-02-10 19:41:00
    推荐学习资料: http://www.cnblogs.com/Ash-ly/p/5452580.html http://ylroki.blog.163.com/blog/static/162978871201032775322518/ ... 一、定义 通过图G的每个节点一次,且...
  • 哈密顿

    2013-11-14 22:21:00
    一个哈密顿环是一条沿着图G的n条边环行的路径,它访问每个结点一次并且返回到它的开始位置。换言之,如果一个哈密顿环在某个结点v1∈V处开始,且G中结点按照v1,v2,…,Vn+l的次序被访问,则边(Vi,Vi+1),1≤i≤n,均...
  • 哈密顿回路及解法

    万次阅读 2017-12-24 15:11:59
    也可以 定义为n+1个相邻顶点v0, v1, … ,vn, v0的一个序列,其中序列的第一个顶点和最后一个顶点是相同的,而其他n-1个顶点是互不相同的。 2、当这个图是加权图时,求该图的最短哈密顿回路,就是传说中的旅行商问题...
  • 哈密顿力学摘要

    2020-12-26 23:08:05
    文章目录前言哈密顿原理经典变分问题数学基础区别函数和泛函变分—自变量不变条件下函数自身的变化变分运算规则积分形式的泛函及其极值问题最速降线问题变分问题的欧拉方程哈密顿原理的优点哈密顿正则方程广义动量与...
  • 矩阵论笔记】零化多项式

    万次阅读 2020-05-07 15:39:42
    矩阵的特征多项式就是矩阵的零化多项式:卡雷哈密顿定理。 特征多项式例: 写成余项的形式,余项要比特征多项式的最高次幂少1阶(不一定1阶) 可以用待定系数求导法来求解余项的系数,求导可以增加式子。要对2阶...
  • 哈密顿路问题

    万次阅读 2013-07-01 21:39:45
    对于一个图中是否存在一条哈密顿路,没有可靠的充分必要条件(貌似邻接矩阵恒式可以?),因此求哈密顿路是一个NP问题,一般要使用搜索和状压dp求解,但汉密尔顿回路的存在有许多充分条件,即当图满足某些特定性质的...
  • \)为哈密顿矩阵\(H\)的本征矢,则有: \[H\left|\psi(t)\right>=E\left|\psi(t)\right> \]此处\(E\)为哈密顿矩阵\(H\)的本征能量,或称为本征值。 绝热演化与量子退火 绝热演化过程可以这...
  • 国内许多工科教材在讲到有关拉普拉斯算子(Δ\DeltaΔ)与哈密顿算子(∇\nabla∇)的内容时含混不清,忽略了许多重要定义,使得一些进一步的推导难以理解 现记录我发现的两个主要问题,并予以解答,希望可以帮助到学习...
  • 我们以非线性薛定ding(NLS)层次结构为例,定义并... 我们从任一Lax矩阵的单峰矩阵构造哈密顿量的显式对偶层级以及触发动力学的Lax表示。 简要介绍了通过连续使用对偶泊松结构来构建Lax对的多维晶格的吸引人的过程。
  • 矩阵 

    千次阅读 2009-06-02 09:30:00
    这个定义很好地解释了Matrix代码制造世界的数学逻辑基础。 数学上,矩阵就是由方程组的系数及常数所构成的方阵。把用在解线性方程组上既方便,又直观。例如对于方程组。 a1x+b1y+c1z=d1 a2x+b2y+c2z=d2 a3x+b3y+...
  • 矩阵

    2008-10-28 21:48:04
    这个定义很好地解释了Matrix代码制造世界的数学逻辑基础。  数学上,矩阵就是由方程组的系数及常数所构成的方阵。把用在解线性方程组上既方便,又直观。例如对于方程组。  a1x+b1y+c1z=d1  a2x+b2y+c2z=d2 ...
  • 前言 时隔一个月,再带大家尝试借助哈密顿环来...先来简单介绍一下哈密顿环的定义(引自维基百科): 哈密顿图是一个无向图,由哈密顿爵士提出, 由指定的起点前往指定的终点,途中经过所有其他节点且只经过一次。 在
  • poj2438 哈密顿

    2016-03-18 19:54:48
    1、哈密顿定义:从一个点出发、经过所有的点必须且只能一次,最终回到起点。图中有边可以不经过,但是不能经过两次。 2、存在的充分条件:图有n个顶点,如果图中任意两个不同的顶点的度数之和大于等于n,则图是...
  • 最近接触了一点雅克比的东西,以前学习雅克比矩阵和雅克比行列式是在... 首先介绍定义,雅克比矩阵是一阶偏导数以一定的方式排列成的矩阵,当其实方阵时,行列式称为雅克比行列式。设有m个n元函数组成的函数组:,称之
  • 先来简单介绍一下哈密顿环的定义(引自维基百科): 哈密顿图是一个无向图,由哈密顿爵士提出, 由指定的起点前往指定的终点,途中经过所有其他节点且只经过一次。 在图论中是指含有哈密顿回路的图, 闭
  • 在之前的一篇文章中,我们介绍过如何将矩阵&概率画成图,读者表示妙不可言。最近,该文作者又动手实现了新的想法:将矩阵画成张量网络图。选自math3ma,作者:Algebra,机器之心编译,参与:李志伟、张倩。今天...
  • 浅谈雅可比矩阵

    万次阅读 多人点赞 2019-04-13 19:30:38
    浅谈雅可比矩阵历史渊源 历史渊源 首先,要先介绍一下——多产堪比欧拉,被广泛认为是历史上三大最具运算能力的数学家之一的雅可比先生。 卡尔·雅可比 1804年12月10日,卡尔·雅可比(Carl Gustav Jacob Jacobi)...
  • 数学 矩阵

    2012-04-23 13:38:00
    数学上,一个m×n的矩阵是一个由m行n列元素排列成的矩形阵列。矩阵里的元素可以是数字、符号或数学式。以下是一个由6个数字符素构成的2行3列的矩阵: 大小相同(行数列数都相同)的矩阵之间可以相互加减,具体是...
  • 相关文件 关注小编,私信小编可以领取源码哟~~ 当然别忘了一键三连哟!...先来简单介绍一下哈密顿环的定义(引自维基百科): 哈密顿图是一个无向图,由哈密顿爵士提出, 由指定的起点前往指定的终点,途中经

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 868
精华内容 347
关键字:

哈密顿矩阵定义