自洽可分的哈密顿系统的辛算法

news2024/11/28 12:51:02
  • 本文只介绍哈密顿系统的辛算法的显式结构
    • 不给出具体的推导过程

自洽可分的哈密顿系统的辛算法

  • 一阶显式辛结构

\begin{matrix} q^{n+1}=q^n+\tau U_p(p^n) \\ p^{n+1}=p^n-\tau V_q(q^{n+1}) \end{matrix}

  • 二阶显式辛结构

\begin{matrix} x=q^n\\ y=p^n-\frac{\tau}{2}V_q(x)\\ q^{n+1}=x+\tau U_q(y)\\ p^{n+1}=y-\frac{\tau}{2}V_q(q^{n+1}) \end{matrix}

  • 四阶显式辛结构

c_1=0,c_2=(2-2^{1/3})^{-1},c_3=1-2(2-2^{1/3}),c_4=(2-2^{1/3})^{-1}

d_1=d_4=(2-2^{1/3})^{-1},d_2=d_3=(1-(2-2^{1/3})^{-1})/2

\begin{matrix} x_1=q^n+c_1\tau U_p(p^n) &y_1 = p^n -d_1V_q(x_1) \\ x_2=x_1+c_2\tau U_p(y_1) & y_2 =y_1-d_2\tau V_q(x_2)\\ x_3 = x_2 + c_3\tau U_p(y_2) &y_3=y_2 -d_3\tau V_q(x_3) \\ q^{n+1}=x_3 +c_4\tau U_p(y_3) &p^{n+1}=y_3-d_4\tau V_q(q^{n+1}) \end{matrix}

  • 全代码
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve

##SymplecticHamilton
##self-consistent
##separable
class symhamcssolver():

    def __init__(self, Up, Vq, 
                 q0=np.array([0]),
                 p0=np.array([1]), 
                 T_start=0, T_end=20, tau=0.01):


        self.Up = Up
        self.Vq = Vq

        self.tau = tau
        self.T = np.arange(T_start, T_end, self.tau)
        self.size = len(self.T)

        self.q = np.zeros((len(self.T), len(q0)))
        self.q[0, :] = q0

        self.p = np.zeros((len(self.T), len(p0)))
        self.p[0, :] = p0
        
        self.tol = 1e-6

    def __str__(self):
        return f"Symplectic algrithm of self-consistent and separable Hamilton system"


    def E1(self,):
        for i in range(1, self.size):
            self.q[i, :] = self.q[i-1, :] + self.tau * self.Up(self.p[i-1, :])
            self.p[i, :] = self.p[i-1, :] - self.tau * self.Vq(self.q[i, :])
        return self.p, self.q
    
    def E2(self,):
        for i in range(1, self.size):
            x = self.q[i-1, :]
            y = self.p[i-1, :] - self.tau/2 * self.Vq(x)
            self.q[i, :] = x + self.tau * self.Up(y)
            self.p[i, :] = y - self.tau/2 * self.Vq(self.q[i, :])
        return self.p, self.q

    def E4(self,):
        alpha = 1/(2-2**(1/3))
        beta  = 1 - 2 * alpha
        c = np.zeros(5)
        d = np.zeros(5)
        c[2] = alpha
        c[4] = alpha
        c[3] = beta
        d[1] = alpha/2
        d[4] = alpha/2
        d[2] = (alpha + beta)/2
        d[3] = (alpha + beta)/2

        for i in range(1, self.size):
            x1 = self.q[i-1, :] + c[1] * self.tau * self.Up(self.p[i-1, :])
            y1 = self.p[i-1, :] - d[1] * self.tau * self.Vq(x1)
            x2 = x1 + c[2] * self.tau * self.Up(y1)
            y2 = y1 - d[2] * self.tau * self.Vq(x2)
            x3 = x2 + c[3] * self.tau * self.Up(y2)
            y3 = y2 - d[3] * self.tau * self.Vq(x3)
            self.q[i, :] = x3 + c[4] * self.tau * self.Up(y3)
            self.p[i, :] = y3 - d[4] * self.tau * self.Vq(self.q[i, :])

        return self.p, self.q


## H = p^2/2m + k/2q^2
m = 1
k = 1
Up = lambda p:p[:]/m
Vq = lambda q:k*q[:]
    
c = symhamcssolver(Up, Vq)

fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(1,3,1)
ax2 = fig.add_subplot(1,3,2)
ax3 = fig.add_subplot(1,3,3)
 
p, q = c.E1()
ax1.plot(p, q)
ax1.set_title("E1")
ax1.set_xlabel("q")
ax1.set_ylabel("p", rotation=0)

p, q = c.E2()
ax2.plot(p, q)
ax2.set_title("E2")
ax2.set_xlabel("q")
ax2.set_ylabel("p", rotation=0) 

p, q = c.E4()
ax3.plot(p, q)
ax3.set_title("E4")
ax3.set_xlabel("q")
ax3.set_ylabel("p", rotation=0)

plt.pause(0.01)




  • 解示意图 

 

 

含时可分哈密顿系统

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
 
##SymplecticHamilton
##self-inconsistent
##separable

##
## \frac{dq}{dt} =  U_p(p, t)
## \frac{dp}{dt} = -V_q(q, t)

##
class symhamicssolver():
 
    def __init__(self, Up, Vq, 
                 q0=np.array([0]),
                 p0=np.array([1]), 
                 T_start=0, T_end=20, tau=0.01):
 
        self.Up = Up
        self.Vq = Vq
 
        self.tau = tau
        self.T = np.arange(T_start, T_end, self.tau)
        self.size = len(self.T)
 
        self.q = np.zeros((len(self.T), len(q0)))
        self.q[0, :] = q0
 
        self.p = np.zeros((len(self.T), len(p0)))
        self.p[0, :] = p0
        
        self.tol = 1e-6
 
    def __str__(self):
        return f"Symplectic algrithm of self-inconsistent and separable Hamilton system"
 
    def E1(self,):
        for i in range(1, self.size):
            self.q[i, :] = self.q[i-1, :] + self.tau * self.Up(self.p[i-1, :], self.T[i-1])
            self.p[i, :] = self.p[i-1, :] - self.tau * self.Vq(self.q[i, :]  , self.T[i])
        return self.p, self.q
    
    def E2(self,):
        for i in range(1, self.size):
            x = self.q[i-1, :]
            y = self.p[i-1, :] - self.tau/2 * self.Vq(x, self.T[i-1])
            self.q[i, :] = x + self.tau * self.Up(y, self.T[i-1]+self.tau/2)
            self.p[i, :] = y - self.tau/2 * self.Vq(self.q[i, :], self.T[i])
        return self.p, self.q
 
    def E4(self,):
        alpha = 1/(2-2**(1/3))
        beta  = 1 - 2 * alpha
        c = np.zeros(5)
        d = np.zeros(5)
        c[2] = alpha
        c[4] = alpha
        c[3] = beta
        d[1] = alpha/2
        d[4] = alpha/2
        d[2] = (alpha + beta)/2
        d[3] = (alpha + beta)/2
 
        for i in range(1, self.size):
            x1 = self.q[i-1, :] + c[1] * self.tau * self.Up(self.p[i-1, :], self.T[i])
            zeta1 = self.T[i-1] + c[1] * self.tau
            y1 = self.p[i-1, :] - d[1] * self.tau * self.Vq(x1, zeta1)
            eta1  = self.T[i-1] + d[1] * self.tau
            
            x2 = x1 + c[2] * self.tau * self.Up(y1, eta1)
            zeta2 = zeta1 + c[2] * self.tau
            y2 = y1 - d[2] * self.tau * self.Vq(x2, zeta2)
            eta2  = eta1 + d[2] * self.tau             
            
            x3 = x2 + c[3] * self.tau * self.Up(y2, eta2)
            zeta3 = zeta2 + c[3] * self.tau
            y3 = y2 - d[3] * self.tau * self.Vq(x3, zeta3)
            eta3  = eta2 + d[3] * self.tau
            
            self.q[i, :] = x3 + c[4] * self.tau * self.Up(y3, eta3)
            self.p[i, :] = y3 - d[4] * self.tau * self.Vq(self.q[i, :], self.T[i])
 
        return self.p, self.q
 
 
## H = p^2/2m + k/2q^2
m = 1
k = 1
Up = lambda p, t:p[:]/m + t
Vq = lambda q, t:k*q[:] + t
    
c = symhamicssolver(Up, Vq)
 
fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(1,3,1)
ax2 = fig.add_subplot(1,3,2)
ax3 = fig.add_subplot(1,3,3)
 
p, q = c.E1()
ax1.plot(p, q)
ax1.set_title("E1")
ax1.set_xlabel("q")
ax1.set_ylabel("p", rotation=0)
 
p, q = c.E2()
ax2.plot(p, q)
ax2.set_title("E2")
ax2.set_xlabel("q")
ax2.set_ylabel("p", rotation=0) 
 
p, q = c.E4()
ax3.plot(p, q)
ax3.set_title("E4")
ax3.set_xlabel("q")
ax3.set_ylabel("p", rotation=0)
 
plt.pause(0.01)
 
 
 
 

 

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1041967.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

Room Arranger for Mac: 轻松创造梦想家园的必备设计软件

你是否曾经梦想过自己动手设计理想中的家居环境?你是否希望通过一个简单易用的工具来实现你的设计理念?那么,Room Arranger for Mac就是你的最佳选择! Room Arranger是一款专门为Mac用户打造的室内设计软件,它拥有直观…

软件测试面试题 —— 整理与解析(4)

😏作者简介:博主是一位测试管理者,同时也是一名对外企业兼职讲师。 📡主页地址:【Austin_zhai】 🙆目的与景愿:旨在于能帮助更多的测试行业人员提升软硬技能,分享行业相关最新信息。…

搭建SpringBoot项目三种方式(超详细版)

目录 一、官网下载压缩包解压 二、通过Idea脚手架搭建 三、Spring Boot项目结构 3.1 pom.xml文件 3.2 启动类 3.3 配置文件 四、通过创建Maven项目添加依赖 一、官网下载压缩包解压 接下来我们搭建一个SpringBoot项目,并引入SpringMVC的功能,首先…

【JAVA EE】详解单点登录

作者简介 目录 1.概述 2.实现方案 2.1.分散鉴权 2.2.集中鉴权 1.概述 SSO,即进行一次认证,然后就可以访问所有子系统。很明显SSO只是一种具象化的目标而已,目前业内为了实现单点登录、统一鉴权,提出了一系列的打法。比如直接…

C# 读取Execl文件3种方法

方法 1,使用OLEDB可以对excel文件进行读取 1.1C#提供的数据连接有哪些 对于不同的.net数据提供者,ADO.NET采用不同的Connection对象连接数据库。这些Connection对我们屏蔽了具体的实现细节,并提供了一种统一的实现方法。 Connection类有四…

ElementUI之首页导航及左侧菜单(模拟实现)

目录 ​编辑 前言 一、mockjs简介 1. 什么是mockjs 2. mockjs的用途 3. 运用mockjs的优势 二、安装与配置mockjs 1. 安装mockjs 2. 引入mockjs 2.1 dev.env.js 2.2 prod.env.js 2.3 main.js 三、mockjs的使用 1. 将资源中的mock文件夹复制到src目录下 2. 点击登…

Java SimpleDateFormat格式化日期时间

java.text.SimpleDateFormat 格式化日期时间, 参考 api 说明 Overview (Java Platform SE 8 ) Examples The following examples show how date and time patterns are interpreted in the U.S. locale. The given date and time are 2001-07-04 12:08:56 local t…

Normalization总结(BN/LN/WN/IN/GN)

一、简介 在深度学习领域,Normalization用得很多,BN(Batch Normalization)于2015年由 Google 提出,开创了Normalization 先河;2016年出了LN(layer normalization)和IN(I…

基于UDP协议的网络服务器的模拟实现

目录 服务端类UdpServer的模拟实现 服务端类UdpServer的成员变量 服务端类UdpServer的构造函数、初始化函数initServer、析构函数 服务端类UdpServer的start函数 服务端类UdpServer的整体代码(即udp_server.h文件的整体代码) 基于服务端类UdpServe…

不同的jdk版本编译得到的class文件中的信息是不是会不一样

不同的jdk版本编译得到的class文件中的信息是不是会不一样 不同的 JDK 版本编译得到的 .class 文件中的信息可能会有所不同。主要的差异可能出现在以下几个方面: 类文件版本号:随着 JDK 版本的升级,类文件的版本号也会发生变化。例如&#x…

左神高级进阶班6(利用快排的partition过程、BFPRT、动态规划的斜率优化技巧、二叉树的递归套路、完美洗牌问题)

目录 【案例1 利用快排的partition过程,BFPRT】 【题目描述】 【思路解析】 【代码实现】 【案例2 动态规划的斜率优化技巧】 【题目描述】 【思路解析】 【代码实现】 【案例3 二叉树的递归套路】 【题目描述】 【搜索二叉树定义】 【思路解析】 【代…

BERT 快速理解——思路简单描述

定义: BERT(Bidirectional Encoder Representations from Transformers)是一种预训练的语言模型,它基于Transformer架构,通过在大规模的未标记文本上进行训练来学习通用的语言表示。 输入 在BERT中,输入…

一篇博客学会系列(1) —— C语言中所有字符串函数以及内存函数的使用和注意事项

目录 1、求字符串长度函数 1.1、strlen 2、字符串拷贝(cpy)、拼接(cat)、比较(cmp)函数 2.1、长度不受限制的字符串函数 2.1.1、strcpy 2.1.2、strcat 2.1.3、strcmp 2.2、长度受限制的字符串函数 2.2.1、strncpy 2.2.2、strncat 2.2.3、strncmp 3、字符串查找函数…

Java 大厂八股文面试专题-JVM相关面试题 垃圾回收算法 GC JVM调优

Java 大厂八股文面试专题-JVM相关面试题 类加载器_软工菜鸡的博客-CSDN博客 3 垃圾收回 3.1 简述Java垃圾回收机制?(GC是什么?为什么要GC) 难易程度:☆☆☆ 出现频率:☆☆☆ 为了让程序员更专注于代码的实现…

如何使用iPhone15在办公室观看家里电脑上的4k电影,实现公网访问本地群晖!

如何使用iPhone15在办公室观看家里电脑上的4k电影? 文章目录 如何使用iPhone15在办公室观看家里电脑上的4k电影?1.使用环境要求:2.下载群晖videostation:3.公网访问本地群晖videostation中的电影:4.公网条件下使用电脑…

【Java基础-JDK21新特性】它发任它发,我用java8

💝💝💝欢迎来到我的博客,很高兴能够在这里和您见面!希望您在这里可以感受到一份轻松愉快的氛围,不仅可以获得有趣的内容和知识,也可以畅所欲言、分享您的想法和见解。 推荐:kuan 的首页,持续学…

黑马JVM总结(二十四)

(1)练习-分析a a:先执行iload1:把数据读入到操作数栈中 iinc:把局部变量表中的1号曹位做一个自增,他在局部变量表中发生的并没有影响到操作数栈 a:限制性自增在做iload 自增变成12 iload把12读取到操作数…

Mac电脑信息大纲记录软件 OmniOutliner 5 Pro for Mac中文

OmniOutliner 5 Pro是一款专业级的Mac大纲制作工具,它可以帮助用户更好地组织和管理信息,以及制作精美的大纲。以下是OmniOutliner 5 Pro的主要功能和特点: 强大的大纲组织和管理功能。OmniOutliner 5 Pro为用户提供了多层次的大纲结构&…

Python语法之条件语句(很详细)

目录 Python条件语句的介绍 定义 if的语法和实例(最基本的) 语法 gif动态图展示 具体实例 实现思路: if-elif-else的语法和实例(最基本的) 语法 具体实例 实现思路: 判断需要多个条件需同时判断语法和实例(最基…

利用Axure RP和cpolar内网穿透实现公网访问本地web网页

AxureRP制作静态站点发布互联网,内网穿透实现公网访问 文章目录 AxureRP制作静态站点发布互联网,内网穿透实现公网访问前言1.在AxureRP中生成HTML文件2.配置IIS服务3.添加防火墙安全策略4.使用cpolar内网穿透实现公网访问4.1 登录cpolar web ui管理界面4…