计算数学
精细积分法
本节介绍由我国钟万勰院士提出的齐次线性定常微分系统的精细积分方法。 精细积分法看上去简单，但是对于这一类微分方程可以获得非常高精度的计算结果，亦既“计算机上的精确解”。 对于不满足精细积分法的方程，常常可以在此基础上以摄动法等手段进行处理。 精细积分法在控制论与弹性力学中已经取得了很重要的应用。下面介绍这一方法。 对于定常齐次微分方程 \[\frac{{d{\mathbf{x}}}}{{dt}} = {\mathbf{Ax}}\] 其通解可以写为 \[{\mathbf{x}}\left( t \right) = \exp \left( {{\mathbf{A}}t} \right){{\mathbf{x}}_0}\] 这里出现了矩阵指数函数，这一概念在大学常微分方程以及矩阵论等课程中都会讲到。 \[\exp \left( {{\mathbf{A}}t} \right) = {{\mathbf{I}}_n} + {\mathbf{A}}t + \frac{{{{\left( {{\mathbf{A}}t} \right)}^2}}}{{2!}} + \frac{{{{\left( {{\mathbf{A}}t} \right)}^3}}}{{3!}} + \cdots \] 设积分步长为\(\eta \)，则 \[{{\mathbf{x}}_1} = {\mathbf{T}}{{\mathbf{x}}_0} = \exp \left( {{\mathbf{A}}t} \right){{\mathbf{x}}_0}\] 由此可以逐步递推 \[{{\mathbf{x}}_1} = {\mathbf{T}}{{\mathbf{x}}_0},{{\mathbf{x}}_2} = {\mathbf{T}}{{\mathbf{x}}_2}, \cdots ,{{\mathbf{x}}_{k + 1}} = {\mathbf{T}}{{\mathbf{x}}_k}, \cdots \] 问题便转到计算矩阵\({\mathbf{T}}\)上来，要求其精确性较高。 矩阵指数计算有两个要点：（1）利用指数函数的加法定理，（2）将注意力集中在增量上，而不是全量。 \[\exp \left( {{\mathbf{A}}\eta } \right) = {\left( {\exp \left( {{\mathbf{A}}\frac{\eta }{m}} \right)} \right)^m}\] 其中$m$为任意整数，这里可以选$$m = {2^N}$$。假设N=20,则$m=1048576$。本来每段区间就不大，则$\tau = \frac{\eta }{m}$ 是一个非常小的区段。在该很小的区段上，幂级数前几项展开已经足够，此时矩阵指数与单位矩阵非常接近，可以写为 \[\exp \left( {{\mathbf{A}}\tau } \right) \approx {{\mathbf{I}}_n} + {{\mathbf{T}}_a}\] \[{{\mathbf{T}}_a} = {\mathbf{A}}\tau + \frac{{{{\left( {{\mathbf{A}}\tau } \right)}^2}}}{2}\left[ {{{\mathbf{I}}_n} + \frac{{{\mathbf{A}}\tau }}{3} + \frac{{{{\left( {{\mathbf{A}}\tau } \right)}^2}}}{{12}}} \right]\] 在计算机中直观重要的一点是矩阵指数的存储只能存储\({{\mathbf{T}}_a}\)，而不是全量。 因为${{\mathbf{T}}_a}$很小，当其与单位矩阵相加时候，其精度将丧失殆尽。这是计算数学中的“大数吃小数现象”，需要特别注意的。 对矩阵做分解 \[{\mathbf{T}} = {\left( {{{\mathbf{I}}_n} + {{\mathbf{T}}_a}} \right)^{2N}} = {\left( {{{\mathbf{I}}_n} + {{\mathbf{T}}_a}} \right)^{2\left( {N - 1} \right)}} \times {\left( {{{\mathbf{I}}_n} + {{\mathbf{T}}_a}} \right)^{2\left( {N - 1} \right)}}\] 这种分解一直做下去，共N次。 在矩阵计算中\({{\mathbf{I}}_b}\)，\({{\mathbf{I}}_c}\)，有 \[\left( {{\mathbf{I}} + {{\mathbf{T}}_b}} \right) \times \left( {{\mathbf{I}} + {{\mathbf{T}}_c}} \right) = {\mathbf{I}} + {{\mathbf{T}}_b} + {{\mathbf{T}}_c} + {{\mathbf{T}}_b} \times {{\mathbf{T}}_c}\] 当\({{\mathbf{I}}_b}\)，\({{\mathbf{I}}_c}\)很小时，不应该加上单位矩阵后再执行乘法，而是将\({{\mathbf{I}}_b}\)，\({{\mathbf{I}}_c}\)都看成\({{\mathbf{I}}_a}\)。 矩阵分解相当于 \[\begin{gathered} {\text{for}}\quad i = 1:N \hfill \\ \quad \quad {{\mathbf{T}}_a} = 2{{\mathbf{T}}_a} + {{\mathbf{T}}_a} \times {{\mathbf{T}}_a} \hfill \\ {\text{end}} \hfill \\ \end{gathered} \] 在做完矩阵N次分解之后，再执行 \[{\mathbf{T}} = {{\mathbf{I}}_b} + {{\mathbf{T}}_a}\] 此时\({{\mathbf{T}}_a}\)已经不是小量，没有严重的舍入误差了，以上便是精细积分法的基本思想与完整过程。下面给出对线性定常系统精细积分方法的MATLAB代码
function T=pre_int(t0,tt,A,x0,step,N)
%---------------------------- General M-file comment
% Version : V1.0
% Coded by : song.yz
% Date : 2012-10-16 13:28:39 *星期二*
%-----------------------------------------------------
% Purpose : 精细积分法
% Reference :
% *. 《应用力学对偶体系》
% *. 《状态空间控制理论与计算》
% *. 《应用力学的辛数学方法》
%-----------------------------------------------------
% Parameters :
% *. t0,tt 积分区间
% *. A-----系数矩阵
% *. x0----------初值
% *. K 区间分为多少，亦为输出步长
% *. N 每段上的划分
% *. T-------T矩阵
%----------------------------------------------------
T=A*2;
fid=fopen('fout_pre_int.txt','w+');
%打开文件
Imax=1e8;
%最大允许积分步数
E=eye(size(A));
M=2^N;
tao=step/M;
At=A*tao;
Ta=At+At*At*(E+At/3+At*At/12)/2;
for i=1:N
Ta=Ta*2+Ta*Ta;
end
T=E+Ta;
%获得T矩阵
x1=x0;
for j=1:Imax
t=t0+j*step;
x2=T*x1;
TMP=[t,x2'];
fprintf(fid,'%10.6f %10.6f %10.6f \n',TMP);
if t>tt
break;
end
x1=x2;
end
status= fclose(fid);
%关闭文件，如果status==0，则关闭成功
%如果status==-1则关闭失败
end