计算数学
常微分方程数值解法
保辛积分
力学系统的哈密顿方程可以通过对力学系统的拉格朗日函数进行勒让德变换得到， 更简便导出正则方程的方法是由相空间的变分原理直接得到。积分曲线$$\int {{\bf{p}}d{\bf{q}} - Hdt} $$当积分路路径端点恒在n维子空间 \(\left( {t = {t_0},{\bf{q = }}{{\bf{q}}_0}} \right)\)上变动时，必定以\(\gamma \)为驻定曲线。 曲线\(\gamma \)是形式\(\left( {t = {t_0},{\bf{q = }}{{\bf{q}}_0}} \right)\)的一条涡线。因此\(\left( {t = {t_0},{\bf{q = }}{{\bf{q}}_0}} \right) \) 在“经过涡方向的无穷小平行四边形”上积分为零。
n个自由度力学系统的哈密顿方程如下 \[\left\{ \begin{array}{l} \frac{{d{\bf{q}}}}{{dt}} = \frac{{\partial H}}{{\partial {\bf{p}}}}\\ \frac{{d{\bf{p}}}}{{dt}} = - \frac{{\partial H}}{{\partial {\bf{q}}}} \end{array} \right.\] 引入2n* 2n矩阵 \[{\bf{J}} = \left[ {\begin{array}{*{20}{c}} {\bf{0}}&{{{\bf{I}}_n}}\\ { - {{\bf{I}}_n}}&{\bf{0}} \end{array}} \right]\] 则哈密顿方程可以写为 \[{\bf{v}} = {\bf{J}}\left( {\frac{{\partial H}}{{\partial {\bf{v}}}}} \right)\] 其中 \[{\bf{v}} = {\left( {{{\bf{q}}^T}{\bf{,}}{{\bf{p}}^T}} \right)^T}\] 在状态向量的形式下，定常正则变换可以用2n维向量的变换描述 \[\varsigma = \varsigma \left( {\bf{v}} \right),{\bf{\varsigma }} = {\bf{\varsigma }}{\left( {{{\bf{Q}}^T},{{\bf{P}}^T}} \right)^T}\] 其逆变换也是定常的，把逆变换和变换分别带入正则方程得到 \[{\bf{\dot \varsigma = }}{\left( {\frac{{\partial {\bf{\varsigma }}}}{{\partial {\bf{v}}}}} \right)^T}{\bf{J}}\left( {\frac{{\partial {\bf{\varsigma }}}}{{\partial {\bf{v}}}}} \right)\frac{{\partial K}}{{\partial {\bf{v}}}}\] 定义 \[{\bf{S}} = {\left( {\frac{{\partial {\bf{\varsigma }}}}{{\partial {\bf{v}}}}} \right)^T}\] 正则变换要求变换后仍然保持正则形式，于是有正则变换的条件 \[{{\bf{S}}^T}{\bf{JS}} = {\bf{J}}\] 满足这一条件的映射称为辛映射，S称为辛矩阵。 辛映射（或辛矩阵）有如下性质：
[封闭性] 辛映射的复合依然是辛映射。即a与b是任意的辛映射， 则有唯一确定的辛映射$$a \circ b = c$$也是辛映射。
[满足结合律] 即对于辛映射a,b,c，满足$$\left( {a \circ b} \right) \circ c = a \circ \left( {b \circ c} \right)$$。
[有单位元素存在] 存在单位映射1，有$$1 \circ a = a = a \circ 1$$。实际上即不进行正则变换，或者说变换就是其本身。
[辛映射存在逆映射] 即对于映射f，总存在逆映射\({{a^{ - 1}}}\)使得$${a^{ - 1}} \circ a = 1 = a \circ {a^{ - 1}}$$。
因此辛映射（或辛矩阵）构成一个群，且有子群，如行列式为1的辛矩阵构成子群。 偶次微分流形的辛结构（symplectic struture中文有人译为辛构造）是一个闭的非退化的微分2形式 \({\omega ^2}\)，在流形$${M^{2n}} = \left\{ {{{\left( {{{\bf{p}}^T}{\bf{,}}{{\bf{q}}^T}} \right)}^T}} \right\}$$中， \({\omega ^2}\)可以表示为 \[{\omega ^2} = d{\bf{p}} \wedge d{\bf{q}}\] 在辛流形\(\left( {{M^{2n}},\omega } \right)\)上，哈密顿函数H的矢量场 \({\bf{J}}dH\)给出一个单参数微分同胚群$${g^t}:{M^{2n}} \to {M^{2n}}$$： \[\frac{d}{{dt}}{|_{t = 0}}{g^t}\left[ \begin{array}{l} {\bf{q}}\\ {\bf{p}} \end{array} \right] = {\bf{J}}dH\left( {\left[ \begin{array}{l} {\bf{q}}\\ {\bf{p}} \end{array} \right]} \right)\] 群\({g^t}\)称为具有哈密顿函数H的哈密顿相流。 哈密顿相流具保持辛结构特性，即 \[{\left( {{g^t}} \right)^*}{\omega ^2} = {\omega ^2}\] 当\(n=1\)时候$${M^{2n}} = {R^2}$$，就是大家熟悉的相流保面积，亦即刘维尔定理。 线性哈密顿系统和可分哈密顿系统都是特殊的哈密顿系统， 可以从一般哈密顿系统的辛格式得到线性哈密顿系统和可分哈密顿系统的辛格式，也可以 针对线性哈密顿系统和可分哈密顿系统的特殊性寻求更为有效的辛格式。 如果哈密顿函数可以表示为 \[{\bf{H}}\left( {\bf{z}} \right){\bf{ = H}}\left( {{\bf{q,p}}} \right){\bf{ = U}}\left( {\bf{p}} \right){\bf{ + V}}\left( {\bf{p}} \right)\] 则称哈密顿系统是可分的，可分哈密顿系统的正则方程为 \[\left\{ \begin{array}{l} \frac{{d{\bf{q}}}}{{dt}} = \frac{{\partial {\bf{U}}\left( {\bf{p}} \right)}}{{\partial {\bf{p}}}}\\ \frac{{d{\bf{p}}}}{{dt}} = - \frac{{\partial {\bf{V}}\left( {\bf{q}} \right)}}{{\partial {\bf{q}}}} \end{array} \right.\] 这里给出一种由日本学者Yoshida构造的4阶显式辛格式 \[{{\bf{x}}^1}{\bf{ = }}{{\bf{p}}^n} + {c_1}h\left( { - \frac{{\partial {\bf{V}}\left( {{{\bf{q}}^n}} \right)}}{{\partial {\bf{q}}}}} \right),{{\bf{y}}^1}{\bf{ = }}{{\bf{q}}^n}{\bf{ + }}{d_1}h\frac{{\partial {\bf{U}}\left( {{{\bf{x}}^1}} \right)}}{{\partial {\bf{p}}}}\] \[{{\bf{x}}^2}{\bf{ = }}{{\bf{x}}^1} + {c_2}h\left( { - \frac{{\partial {\bf{V}}\left( {{{\bf{y}}^1}} \right)}}{{\partial {\bf{q}}}}} \right),{{\bf{y}}^2}{\bf{ = }}{{\bf{y}}^1}{\bf{ + }}{d_2}h\frac{{\partial {\bf{U}}\left( {{{\bf{x}}^2}} \right)}}{{\partial {\bf{p}}}}\] \[{{\bf{x}}^3}{\bf{ = }}{{\bf{x}}^2} + {c_3}h\left( { - \frac{{\partial {\bf{V}}\left( {{{\bf{y}}^2}} \right)}}{{\partial {\bf{q}}}}} \right),{{\bf{y}}^3}{\bf{ = }}{{\bf{y}}^2}{\bf{ + }}{d_2}h\frac{{\partial {\bf{U}}\left( {{{\bf{x}}^3}} \right)}}{{\partial {\bf{p}}}}\] \[{{\bf{p}}^{n + 1}}{\bf{ = }}{{\bf{x}}^3} + {c_4}h\left( { - \frac{{\partial {\bf{V}}\left( {{{\bf{y}}^3}} \right)}}{{\partial {\bf{q}}}}} \right),{{\bf{q}}^{n + 1}} = {{\bf{y}}^3}{\bf{ + }}{d_4}h\frac{{\partial {\bf{U}}\left( {{{\bf{p}}^{n + 1}}} \right)}}{{\partial {\bf{p}}}}\] 其中系数可以取 \[\left( \begin{array}{l} {c_1}\\ {c_2}\\ {c_3}\\ {c_4} \end{array} \right) = \left( \begin{array}{l} 0\\ \alpha \\ \beta \\ \alpha \end{array} \right),\quad \quad \left( \begin{array}{l} {d_1}\\ {d_2}\\ {d_3}\\ {d_4} \end{array} \right) = \left( \begin{array}{l} \frac{\alpha }{2}\\ \frac{{\alpha + \beta }}{2}\\ \frac{{\alpha + \beta }}{2}\\ \frac{\alpha }{2} \end{array} \right)\] 或者取 \[\left( \begin{array}{l} {c_1}\\ {c_2}\\ {c_3}\\ {c_4} \end{array} \right) = \left( \begin{array}{l} \frac{\alpha }{2}\\ \frac{{\alpha + \beta }}{2}\\ \frac{{\alpha + \beta }}{2}\\ \frac{\alpha }{2} \end{array} \right),\quad \quad \left( \begin{array}{l} {d_1}\\ {d_2}\\ {d_3}\\ {d_4} \end{array} \right) = \left( \begin{array}{l} \alpha \\ \beta \\ \alpha \\ 0 \end{array} \right)\] 这里$$\alpha = {\left( {2 - {2^{\frac{1}{3}}}} \right)^{ - 1}},\beta = 1 - 2\alpha $$。 辛算法在长时间、多步数的计算中和保持系统整体结构上较传统的非辛算法显示出了 强大的优越性。目前，辛算法在天体力学的演化问题，分子动力学和量子物理等问题上取得了较多的研究成果。 本节以一个简单的范例对辛算法与普通的非辛算法进行积分比较，如果欲深入系统的了解 其他的辛差分格式请参考冯康、秦孟兆所著的《哈密顿系统的辛几何算法》。 谐振子是许多学科中最基本的问题，如量子力学、声学与振动等。 对于理想谐振子系统的拉格朗日函数为 \[L = T - V = \frac{1}{2}m\dot x - \frac{1}{2}k{x^2}\] 广义动量 \[p = \frac{{\partial L}}{{\partial \dot x}} = m\dot x\] 哈密顿函数为 \[H = p\dot x - L = \frac{1}{m}{p^2} + \frac{1}{2}k{x^2}\] 由此可以建立哈密顿正则方程 \[\left\{ \begin{array}{l} \frac{{dx}}{{dt}} = \frac{p}{m}\\ \frac{{dp}}{{dt}} = - kx \end{array} \right.\] 这里以谐振子的长期演化来检验辛算法。代码为C#版本，代码用到C#定义的一些数学类，其定义方法见《C#科学计算讲义》。
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
class cHamilton
{
/*--------------------------------------class comment
Version : V1.0
Coded by : syz
Date : 2011-06-28 16:19:05 *星期二*
----------------------------------------------------
parameters :
* k--弹簧劲度系数
* m--谐振子质量
*
Methods :
* Hamilton ----辛格式计算谐振子问题
* RK4----------4阶RK方法解算谐振子
* func---------谐振子的正则方程（哈密顿方程）
--------------------------------------------------*/
public double k = 1.0;
//弹簧劲度系数，无量纲化
public double m = 1.0;
//振子质量，无量纲化
public void Hamilton(long NN, double t0, double p0, double q0, double h)
{
StreamWriter ob1 = new StreamWriter("fout.txt");
// ob1.WriteLine("辛RK方法");
double afa = 1.0 / (2.0 - Math.Pow(2.0, 1.0 / 3));
double beta = 1 - 2.0 * afa;
double c1 = 0.0;
double c2 = afa;
double c3 = beta;
double c4 = afa;
double d1 = afa / 2.0;
double d2 = (afa + beta) / 2.0;
double d3 = d2;
double d4 = d1;
double x1, x2, x3, y1, y2, y3, p, q, p1, q1;
//积分初值
p = p0;
q = q0;
for (int i = 0; i < NN; i++)
{
x1 = p - c1 * h * k*q;
y1 = q + d1 * h * p / m;
x2 = x1 - c2 * h * k * y1;
y2 = y1 + d2 * h * x2 / m;
x3 = x2 - c3 * h * k * y2;
y3 = y2 + d3 * h * x3 / m;
p1 = x3 - c4 * h * k * y3;
q1 = y3 + d4 * h * p1 / m;
p = p1;
q = q1;
//输出广义坐标与广义动量
ob1.WriteLine("{0:F8} {1:F8}", p, q);
}
ob1.Close();
}
public void RK4(int N, int M, double t0, double h, double[] y0)
{
/* ----------------------------------------------------
Desciption :
parameters :
* N----------积分步数
* t0----积分初始时刻
* h-----积分步长
* y0[M]------积分初值
* M----------方程维数
Methods :
* solve----RK4方法函数
* func-----求解的方程函数
*
--------------------------------------------------*/
//生成输出文件
StreamWriter ob2 = new StreamWriter("fRK.txt");
// ob2.WriteLine("RK4方法解微分方程组");
int i;
//积分的第i步
double[] k1 = new double[M];
double[] k2 = new double[M];
double[] k3 = new double[M];
double[] k4 = new double[M];
//注意这里不能再初始化，否则y0就会被重新设置为0
// y0=new double[M];
//初值向量
//计算结果，用文件输出
double[] y = new double[M];
// double h, t;
// h = (tt - t0) / N;
//步长
//积分初值
y = y0;
double t = t0;
double[] vec1 = new double[M];
double[] vec2 = new double[M];
double[] vec3 = new double[M];
int j; //分量指标
for (i = 0; i < N; i++)
{
//第一次调用函数
func(out k1, y, t);
for (j = 0; j < M; j++)
vec1[j] = y[j] + 0.5 * h * k1[j];
//第二次调用函数
func(out k2, vec1, t + 0.5 * h);
for (j = 0; j < M; j++)
vec2[j] = y[j] + 0.5 * h * k2[j];
//第三次调用函数
func(out k3, vec2, t + 0.5 * h);
for (j = 0; j < M; j++)
vec3[j] = y[j] + h * k3[j];
//第四次调用函数
func(out k4, vec3, t + h);
for (j = 0; j < M; j++)
y[j] = y[j] + (k1[j] + k2[j] * 2 + k3[j] * 2 + k4[j]) * h / 6.0;
t = t + h;
//输出计算结果
ob2.WriteLine("{0:F8} {1:F8}", y[0], y[1]);
}
ob2.Close();
//关闭文件
}
public void func(out double[] f, double[] y, double t)
{
//右函数
//y应变量
//t自变量
//f方程函数
//不能对y重新执行以下语句，否则参数传递进来的y将被置于0
//y = new double[2];
f = new double[2];
f[0] = y[1]/m;
f[1] = -k*y[0];
}
}
namespace Hamilton
{
class Program
{
static void Main(string[] args)
{
cHamilton ob;
ob=new cHamilton();
//用辛算法解谐振子问题
ob.Hamilton(100000,0.0, 20.0, 8.0, 0.5);
//Hamilton(long NN, double t0,double p0,double q0, double h)
double[] y0 = new double[2];
y0[0] = 20.0;
y0[1] = 8.0;
ob.RK4(10000, 2, 0.0, 0.5, y0);
//RK4(int N, int M, double t0, double h, double[] y0)
}
}
}