MATHEMATICA 教程

NDSolve 的"OrthogonalProjection"方法

引言

考虑该矩阵微分方程:

其中给出初值 . 假设 ,该解具有保持标准正交性的属性,,并且它对于所有 满秩.

从数值的角度来看,一个关键的问题是如何对一个正交矩阵微分系统数值积分,积分的方式使得数值解仍然正交. 有一些可能的策略. 一种方法,在 [DRV94] 中提出,是采用隐式朗格-库塔方法(如高斯方案). 一些替代的策略在 [DV99] 和 [DL01] 中有描述.         

这里采用的方法是使用任何合理的数值积分方法,然后在每个积分步的末尾使用投影程序进行后处理.

该实现方法的一个重要特点是基本积分方法可以是任何内置的数值计算方法,或者甚至是一个用户定义的过程. 在下面的例子中,一个显式朗格-库塔方法用于基本时间步进(the basic time stepping). 但是,如果需要更高的精确度,则可以方便的使用一个外推法( extrapolation method),例如,通过简单地设置适当的 Method 选项.

投影步(Projection Step)

在每个数值积分步的末尾,用户需要转换该微分系统的近似解矩阵以获取一个正交矩阵. 这可以通过一些方式实现(例如,见 [DRV94] 和 [H97]):

  • 牛顿或者 Schulz 迭代
  • QR 分解
    • 奇异值分解

    牛顿和 Schulz 方法二次收敛,而迭代次数可能会因为在数值积分中采用的误差容差(error tolerances)而有所不同. 在IEEE双精度运算下,一个或者两个迭代对于收敛到正交极性因子通常是足够的(如下所示).

    QR 分解比奇异值分解代价更低,大约是两倍的因子( a factor of two),但是它给出的不是最接近的可能投影.

    定义 (Thin 奇异值分解 [GVL96]):给定一个矩阵 其中 ,存在两个矩阵 使得 的奇异值的对角矩阵,,其中 . 具有标准正交列,而 是正交化的.

    定义 (极分解):给定一个矩阵 和它的奇异值分解 的极分解由两个矩阵 的积给出,其中 并且 . 具有标准正交列,而 是对称半正定的.

    的标准正交极因子 是求解下列式子的矩阵:

    对于 2 和 Frobenius 范数 [H96].

Schulz 迭代

选择的方法基于 Schulz 迭代,它对于 直接适用. 相比之下, 迭代需要先进行 QR 分解.

与基于奇异值分解的直接计算的比较也在这里给出.

Schulz 迭代由下式给出:

Schulz 迭代在每个 浮点运算的迭代中具有一个运算操作,但是具有更丰富的矩阵乘法 [H97].

在一个实际实现中,LAPACK [LAPACK99] 的GEMM-based level 3 BLAS 可以通过自动调谐线性代数软件 [ATLAS00],与特定体系的优化结合起来使用. 这样的考虑意味着,Schulz 迭代的运算操作数目不一定是观察到的计算代价的准确反映. 从 的标准正交性的偏离的有用的界限在 [H89] 中给出:. 与 Schulz 迭代的比较给出对于一些容差 的终止标准 .

标准公式

假定给出常微分方程的当前解的一个初值 ,以及从一步数值积分方法所得的一个解 . 假定控制Schulz 迭代的一个绝对容差 也预先给定.

下面算法可以用于实现.

第 1 步. 设置 .

第 2 步. 计算 .

第 3 步. 计算 .

第 4 步. 如果 或者 , 那么返回 .

第 5 步. 设置 并且跳到第2步.

增量公式

NDSolve 使用赔偿总和来降低由于在每个积分步重复增加小数值 产生的舍入误差的效果 [H96]. 因此,增量 由基本积分器返回.

投影迭代的一个合适的正交修正 可以使用下面算法确定.

第 1 步. 设置 .

第 2 步. 设置 .

第 3 步. 计算 .

第 4 步. 计算 .

第 5 步. 如果 或者 ,那么返回 .

第 6 步. 设置 并且跳到第 2 步.

该修正算法在 中使用,并且显示了使用迭代算法比使用直接过程的优点,由于一个正交修正如何对直接方法推导并不明显.

示例

正交误差测量

计算一个矩阵 的Frobenius 范数 的一个函数可以关于下面的 Norm 函数定义.
In[1]:=
Click for copyable input
的标准正交性偏离的上界可以使用该函数 [H97] 测量.
In[2]:=
Click for copyable input
这里定义功用函数,用于可视化在一个数值积分中的正交误差.
In[4]:=
Click for copyable input
In[6]:=
Click for copyable input

方形系统(Square Systems)

这个例子涉及一个矩阵微分系统在正交群 上的解(见 [Z98]).

矩形微分系统由下列式子给出

其中

并且

解的演化为:

的特征值是 , , . 因此,当 接近 时, 的两个特征值接近 -1. 在区间 上执行数值积分.
In[7]:=
Click for copyable input
In[8]:=
Click for copyable input
这里使用显式朗格-库塔方法计算解. 适当的初始步长和方法的阶数被自动选择,并且在积分区间中步长可能有所不同,步长的选择为了满足局部相对和绝对误差容差. 另一方面,方法的阶数可以使用一个 Method 选项指定.
In[16]:=
Click for copyable input
这里当积分进行时,计算正交误差,或者从正交流形的绝对偏差. 该误差是数值方法的局部准确度的阶数.
In[17]:=
Click for copyable input
Out[18]=
这里使用用于基本积分步骤的显式朗格-库塔法的正交投影法,计算解. 初始步长和方法的阶数与之前的相同,但是在积分中的步长序列可能有所不同.
In[19]:=
Click for copyable input
利用正交投影法,正交误差被降低为接近IEEE双精度运算中四舍五入的水平.
In[20]:=
Click for copyable input
Out[21]=
Schulz 迭代,利用增量公式,一般产生比直接奇异值分解更小的误差.

矩形系统

下面的例子表明正交投影方法的实现对于矩形矩阵微分系统也适用. 正式地表示为,我们感兴趣的是在 Stiefel 流形 上,求解常微分方程,即 × 正交矩阵的集合,其中 .

定义 × 正交矩阵的 Stiefel 流形 是集合 , , 其中 × 单位矩阵.

在 Stiefel 流形上演化的解找到许多应用,如在数值线性代数中的特征值问题,动力系统和信号处理的 Lyapunov 指数的计算.

考虑从 [DL01] 中采用的一个例子:

其中 , ,而 并且 .

精确解由下列式子给出:

标准化得到:

因此, 满足下列 上的弱斜对称系统:

在下面的例子中,该系统在区间 中被求解,其中 并且维度 .
In[22]:=
Click for copyable input
这里计算精确解,该解可以在积分区间内被计算.
In[36]:=
Click for copyable input
这里使用显式朗格-库塔方法计算解.
In[37]:=
Click for copyable input
这里计算在每个积分区间末尾的按分量的绝对全局误差.
In[39]:=
Click for copyable input
Out[39]=
这里计算正交误差——一个 Stiefel 流形偏差的测量.
In[40]:=
Click for copyable input
Out[40]=
这里使用一个正交投影法计算解,该方法使用朗格-库塔方法作为基本的数值积分方案.
In[41]:=
Click for copyable input
在积分区间结尾处的,按分量的绝对全局误差(The componentwise absolute global error)大致与之前的相同,这是因为在数值积分中使用的绝对和相对容差是相同的.
In[43]:=
Click for copyable input
Out[43]=
利用正交投影法,但是,Stiefel 流形的偏差降低到四舍五入的水平.
In[44]:=
Click for copyable input
Out[44]=

实现

方法 的实现具有三个基本组成部分:

  • 初始化. 建立积分中使用的基本方法,决定任何方法系数,并且建立应该使用的任何工作区(workspaces). 在任何实际积分进行前,这只执行一次,并且所得的 对象被验证,以便在每个积分步,不需要被检查. 在这个阶段,对系数的维度和初始条件的一致性进行检查.
  • 在每一步调用基本数值积分方法.
    • 执行正交投影. 这里执行各种测试,如检查基本技法正确进行,并且 Schulz 迭代是收敛的.

    选项可以用来修改 Schulz 迭代的终止标准. 由代码提供的一个选项是 ,它允许对迭代的容差 进行控制. 因子与 Last Place 中的一个单位结合起来,根据积分中采用的工作精度确定(对于IEEE 双精度,有 ).

    终止标准所使用的 Frobenius 范数可以使用 LAPACK LANGE 函数 [LAPACK99] 有效地计算.

    选项 MaxIterations 控制应该执行的最大迭代次数.

选项总结

选项名
默认值
Dimensions{}指定矩阵微分系统的维度
"IterationSafetyFactor"指定 Schulz 迭代 (1) 的终止标准中所采用的安全因子
MaxIterationsAutomatic指定在 Schulz 迭代 (2)中采用的最大迭代次数
Method"StiffnessSwitching"指定数值积分采用的方法

方法 的选项.

New to Mathematica? Find your learning path »
Have a question? Ask support »