精细积分法文献恢复

Posted on
2020-06-30 11:34
hwzhou
阅读(…)
评论(…)
编辑
 收藏

主要文献: 谭述君, 高强, 钟万勰. Duhamel项的精细积分方法在非线性微分方程数值求解中的应用[J]. 计算力学学报, 2010(05):13-19

推导过程

控制方程及其一般解

采用状态空间，即令${\bf p}={x_1,x_2,…,x_L,\dot x_1,\dot x_2,…,\dot x_L}$，则可将振动方程描述为一阶常微分方程组：

\[\left\{ \begin{array}{l}
\bf \dot p=H\cdot p+g(p,t)\\
\bf p(0)=p_0
\end{array} \right.\label{eq:Governing EQ}
\]

其中，$\bf g$为非齐次项，通常是由外力以及非线性恢复力组成。方程$\eqref{eq:Governing EQ} $ 的一般解为：

\[{\bf p}(t)=\exp({\bf H}t)\cdot {\bf p}_0+\int_0^t\exp[{\bf H}\cdot(t-\tau)]\times{\bf g}({\bf p}(\tau),\tau){\rm d}\tau
\]

其中，第一项为通解，齐次项；第二项为Duhamel积分，非齐次项。

将时间域等间隔离散：$t=\{t_1,t_2,…,t_k,…\}$,记步长为$\eta=t_{k+1}-t_k$，$t+t_k$时的函数值为${\bf p}_k$。若${\bf p}_k$已精确求出，则在时间区间$[t_k,t_{k+1}]$内，可写出其递推公式为：

\[{\bf p}_{k+1}=\exp({\bf H}\eta)\cdot {\bf p}_k+\int_0^\eta\exp[{\bf H}\cdot(\eta-\tau)]\cdot{\bf g}({\bf p}(t_k+\tau),t_k+\tau){\rm d}\tau
\label{eq:pim recursion}
\]

齐次项的精细积分法

记${\bf T}(\eta)=\exp({\bf H}\eta)$，为步长$\eta$时，通解（齐次项）的状态传递矩阵。将$T$泰勒展开为：

\[{\bf T}(\eta)={\bf I}+\sum\limits_{j=1}^\infty \frac{{\bf H}^j\cdot \eta^j}{j!}={\bf I}+{\bf T}_a(\eta)
\label{eq:T eta}
\]

用上式计算${\bf T}(\eta)$缺点是,要取足够多的项才能获得所需精度,特别是在$\bf H \eta$的范数$\bf ||H\eta||$较大时,尤为困难。钟教授在“钟万勰.结构动力方程的精细时程积分法.大连理工大学学报,1994;34(2)“中提出的方法,利用了指数函数性质$\exp({\bf H \eta)=\exp (H\eta} /M)^M$,利用尺度因子$m$作用,先计算$\exp[{\bf H} \delta]$，其中$\delta=\eta/M,M=2^N$,再经过$N$次矩阵自乘,解决了上述问题^[1]。

在计算${\bf T}(\delta)$时，由于增量阵${\bf T}_a(\delta)$很小，若在计算机里与单位矩阵直接相加极可能会因大数加小数而损失有效数字。利用其主量不变(单位矩阵$\bf I$)的性质，将增量阵单独存储，并按如下方法计算出步长$\eta$时的增量阵${\bf T}_a(\eta)$（非小量），再将其与主量相加，得到状态传递矩阵的全量${\bf T}(\eta)$。

增量阵加法定理：

\[{\bf T(2\delta)}=\left({\bf I}+{\bf T}_a(\delta)\right)^2={\bf I}+2{\bf T}_a(\delta)+{\bf T}_a(\delta)^2={\bf I}+{\bf T}_a(2\delta)
\]

即：

\[{\bf T}_a(2\delta)=2{\bf T}_a(\delta)+{\bf T}_a(\delta)^2
\label{eq:Ta2delta}
\]

可迭代$N$次，计算出${\bf T}_a(\eta)$如

\[{\bf T}_a(\delta)\rightarrow {\bf T}_a(2^1*\delta) \rightarrow {\bf T}_a(2^2*\delta) \rightarrow…\rightarrow {\bf T}_a(2^N*\delta)={\bf T}_a(\eta)
\label{eq:pim iteration}
\]

非齐次项的精细积分法

外力（含非线性内力）${\bf g(p},t)$在时间区间$[t_k,t_{k+1}]$内采用$\overline m$阶多项式插值

\[{\bf g(p,}t_k+\tau)={\bf g}_{0,k}+{\bf g}_{1,k}\tau+…++{\bf g}_{\overline m,k}\tau^\overline m,\tau\in[0,\eta]
\]

其中，${\bf g}_{m,k},m=1,2,…\overline m$为系数，与时间$\tau$无关。

将上式及式$\eqref{eq:T eta}$代入式$\eqref{eq:pim recursion}$，得：

\[{\bf p}_{k+1}={\bf T}(\eta)\cdot{\bf p}_k+\sum\limits_{m=0}^\overline m \eta^{m+1} \cdot{\bf T}_{m}(\eta)\cdot{\bf g}_{m,k}
\label{eq:pim pk+1}
\]

其中

\[{\bf T}_m(\eta)=\frac{1}{\eta^{m+1}}\int_0^\eta\exp[{\bf H}\cdot(\eta-\tau)]\cdot \tau^m{\rm d}\tau,\hspace{0.5em}m=0,1,2,…,\overline m
\label{eq:Ti eta}
\]

称为：外力（含非线性内力）的响应矩阵或Duhamel积分矩阵。

与齐次项的传递矩阵${\bf T}(\eta)$类似，将${\bf T}_m(\eta)$泰勒展开(原文献中笔误）：

\[\begin{aligned}
{\bf T}_m(\eta)&=\frac{1}{m+1}{\bf I}+\frac{1}{\eta^{m+1}}\sum\limits_{j=1}^\infty\frac{{\bf H}^j}{j!}\int_0^\eta(\eta-\tau)^j\cdot\tau^m{\rm d}\tau\\
&=\frac{1}{m+1}{\bf I}+{\bf T}_{m,a}(\eta)
\end{aligned}
\]

记${\bf T}_{m,a}(\delta)=\sum\limits_{j=1}^\infty {\bf T}_{m,j,a}(\delta)$,则：

\[\begin{equation}\begin{aligned}{\bf T}_{m,j,a}(\delta)&=\frac{{\bf H}^j}{\delta^{m+1}\cdot j!}\int_0^\delta(\delta-\tau)^j\cdot\tau^m{\rm d}\tau\\&=\frac{{\bf H}^j\cdot\delta^j}{\delta^{m+1}\cdot j!}\int_0^\delta \sum\limits_{k=0}^j \C_j^k \delta^k (-\tau)^{j-k}\cdot \tau^m{\rm d}\tau\\&=\frac{{\bf H}^j}{\delta^{m+1}\cdot j!}\sum\limits_{k=0}^j (-1)^{j-k}\delta^k\C_j^k \int_0^\delta\tau^{j-k+m}{\bf d}\tau\\&=\frac{{\bf H}^j}{\delta^{m+1}\cdot j!}\sum\limits_{k=0}^j (-1)^{j-k}\delta^k\C_j^k \frac{\delta^{j-k+m+1}}{j-k+m+1}\\&=\sum\limits_{k=0}^j \frac{(-1)^{j-k}\C_j^k}{(j-k+m+1)\cdot j!}{\bf H}^j\cdot\delta^j\\&=\alpha_{j,m}\hspace{0.2em}{\bf H}^j\delta^j\end{aligned}\end{equation}
\]

而：

\[\begin{equation}\begin{aligned}{\bf T}_{m,a}(\delta)&=\sum\limits_{j=1}^\infty \alpha_{j,m}{\bf H}^j\delta^j\hspace{0.2em},m=0,1,2,…,\overline m\\
\alpha_{j,m}&=\sum\limits_{k=0}^j \frac{(-1)^{j-k}\C_j^k}{(j-k+m+1)\cdot j!}\end{aligned}\end{equation}
\]

其中,$\C_j^k=\dfrac{j!}{k!(j-k)!}$ 为组合数。

类似的，在精细区间$\delta$内的增量${\bf T}_{m,a}(\delta)$依然是小量，需要单独存储，并按类似$\eqref{eq:pim iteration}$的方式迭代至${\bf T}_{m,a}(\eta)$，其关键在于推导出${\bf T}_{m,a}(\delta)$的加法定理（类似式$\eqref{eq:Ta2delta}$，即将${\bf T}_{m,a}(2\delta)$表示成${\bf T}_{m,a}(\delta)$的多项式）。

${\bf T}_{m,a}(\delta)$的加法定理

根据式$\eqref{eq:Ti eta}$，${\bf T}_{m}(2\delta)$可写为：

\[\begin{equation}\begin{aligned}
{\bf T}_m(2\delta)&=\frac{1}{(2\delta)^{m+1}}\int_0^{2\delta}\exp[{\bf H}\cdot(2\delta-\tau)]\cdot \tau^m{\rm d}\tau\\
&=\frac{1}{(2\delta)^{m+1}}\left[\int_0^{\delta}\exp[{\bf H}\cdot(2\delta-\tau)]\cdot \tau^m{\rm d}\tau+\int_\delta^{2\delta}\exp[{\bf H}\cdot(2\delta-\tau)]\cdot \tau^m{\rm d}\tau\right]\\
&=\frac{1}{(2\delta)^{m+1}}\left[\bf A+B\right]
\end{aligned}\end{equation}
\]

其中，$A$和$B$分别计算如下。

\[\begin{equation}\begin{aligned}
{\bf A}&=\frac{1}{(2\delta)^{m+1}}\int_0^{\delta}\exp({\bf H}\delta)\exp[{\bf H}\cdot(\delta-\tau)]\cdot \tau^m{\rm d}\tau\\
&=\frac{\exp({\bf H}\delta)}{2^{m+1}}\cdot \frac{1}{\delta^{m+1}}\int_0^{\delta} \exp[{\bf H}\cdot(\delta-\tau)]\cdot \tau^m{\rm d}\\
&=\frac{1}{2^{m+1}}{\bf T}(\delta)\cdot{\bf T}_i(\delta)\\
&=\frac{1}{2^{m+1}}({\bf I}+{\bf T}_a(\delta))\cdot\left(\frac{\bf I}{m+1}+{\bf T}_{m,a}(\delta)\right)\\
\end{aligned}\end{equation}
\label{eq:Ti2delta A}
\]

\[\begin{equation}\begin{aligned}
{\bf B}&=\frac{1}{(2\delta)^{m+1}}\int_0^{\delta}\exp[{\bf H}\cdot(2\delta-(x+\delta)]\cdot (x+\delta)^m{\rm d}x\\
&=\frac{1}{(2\delta)^{m+1}}\int_0^{\delta}\exp[{\bf H}\cdot(\delta-x)]\cdot \sum\limits_{k=0}^{m}\C_ m ^k\delta^{k} x^{m-k}{\rm d}x\\
&=\sum\limits_{k=0}^{m}\frac{\C_m^k\delta^{k}}{2^{m+1}}\cdot\frac{1}{\delta^{m+1}}\int_0^{\delta}\exp[{\bf H}\cdot(\delta-x)]\cdot x^{m-k}{\rm d}x\\
&=\sum\limits_{k=0}^{m}\frac{\C_ m ^k}{2^{m+1}}\cdot\frac{1}{\delta^{m-k+1}}\int_0^{\delta}\exp[{\bf H}\cdot(\delta-x)]\cdot x^{m-k}{\rm d}x\\
&=\frac{1}{2^{m+1}}\sum\limits_{k=0}^{m}\C_ m ^k\cdot{\bf T}_{m-k}(\delta)\\
&=\frac{1}{2^{m+1}}\sum\limits_{k=0}^{m}\C_m^k\cdot\left(\frac{\bf I}{m-k+1}+{\bf T}_{m-k,a}(\delta)\right)\\
\end{aligned}\end{equation}
\label{eq:Ti2delta B}
\]

则，可分别计算${\bf T}_{i}(2\delta)$（即$\bf A+B$）的主量和增量。

\[\begin{equation}\begin{aligned}
{\bf A+B}|_{主量}&=\frac{1}{2^{m+1}}\left(\frac{1}{m+1}+\sum\limits_{k=0}^m \frac{\C_m^k}{m-k+1}\right){\bf I}\\
&=\frac{1}{2^{m+1}}\left(\frac{1}{m+1}+\sum\limits_{k=0}^m \frac{m!}{(m-k+1)(m-k)!k!}\right){\bf I}\\
&=\frac{1}{2^{m+1}}\left(\frac{1}{m+1}+\sum\limits_{k=0}^m \frac{1}{m+1}\frac{(m+1)!}{(m+1-k)!k!}\right){\bf I}\\
&=\frac{1}{2^{m+1}}\left[\frac{1}{m+1}\left(1+\sum\limits_{k=0}^m \C_{m+1}^k\right)\right]{\bf I}\\
&=\frac{1}{2^{m+1}}\left(\frac{1}{m+1}\cdot2^{m+1}\right){\bf I}\\
&=\frac{1}{m+1}{\bf I}\\
\end{aligned}\end{equation}
\]

\[\begin{equation}\begin{aligned}
{\bf T}_{m,a}(2\delta)&={\bf A+B}|_{增量}\\
&=\frac{1}{2^{m+1}}\left[\frac{{\bf T}_a(\delta)}{m+1}+{\bf T}_{m,a}(\delta)+{\bf T}_{a}(\delta){\bf T}_{m,a}(\delta)+\sum\limits_{k=0}^m\C_m^k{\bf T}_{m-k,a}(\delta) \right]\\
&=\frac{1}{2^{m+1}}\left[\frac{{\bf T}_a(\delta)}{m+1}+2{\bf T}_{m,a}(\delta)+{\bf T}_{a}(\delta){\bf T}_{m,a}(\delta)+\sum\limits_{k=1}^m\C_m^k{\bf T}_{m-k,a}(\delta) \right]
\end{aligned}\end{equation}
\]

迭代过程：

设置增量矩阵数组：${\bf T}_a|_{2L\times 2L}$和${\bf T}_{i,a}|_{2L\times 2L\times i}$，其中$L$为参与计算的模态阶数。

Step 1: $\dfrac{{\bf T}_a}{i+1}+2{\bf T}_{i,a}(:,:,i)+{\bf T}_{a}{\bf T}_{i,a}(:,:,i)\rightarrow {\bf T}_{i,a}(:,:,i)$

Step 2: for $k=i,1,-1$，

一般的，但$N$取10或20时，$\bf ||H\delta||$很小。 ↩︎

本文链接：https://www.cnblogs.com/hwzhou/p/13212794.html

精细积分法文献恢复