求解病态系统的Rosenbrock算法实现
EmilMatthew (EmilMatthew@126.com)
[ 类别 ]算法实现
[推荐指数]★★★
[ 摘要 ]本文主要介绍了一种求解病态系统较为有效的半稳式Range-Kuta算法Rosenbrock的基本原理及其实现。
[ 关键词 ]半隐式Range-Kuta算法，Rosenbrock算法
The Implementation of Rosenbrock Algorithm------
A method for solving stiff system
[Classify] Algorithm Implementation
[ Level ] ★★★
[Abstract] In this article, I mainly introduce the benchmark idea of a Semi Implicit Range-Kuta algorithm, which named as Rosenbrock algorithm, and also with its implementation.
[Key Words]Semi Implicit Range-Kuta, Rosenbrock Algorithm
[0引言]
常微分方程组(Ordinary Differential Equations)的数值解法理论，在工程及数值实现上有着非常重要的应用，其中有一类常微分方程组，其系统呈现病态特性（或称刚性，此类系统常称作stiff system,简介参附录1），对普通的数值求解算法的步长设定要求较高，往往需要设定非常小的步长才能得出满意结果。在一些需要时实计算的领域，如控制系统等，是不允许在求解的时间上有过多耗费的，这就要求使用能在较短时间内就可得出计算结果的数值算法。
[1算法介绍]
由于求解普通常微分方程组的Range-Kuta算法具有求解精度高，适用范围广的特点，因此，针对病态系统的数值算法构造有许多也是基于Range-Kuta算法展开的。计算数学工作者发现，如将Range-Kuta算法表达成隐式的形式（参附录2），并对其进行求解，则可以在较短的步数内求得病态系统的状态变量值。
但是，如果直接对隐式Range-Kuta方程进行求解，将面临着耗时、迭代易发散的不良特性。因此，真正有效的求解病态系统的数值算法，有很大一部分基于对隐式Range-Kuta方程展开的。这里要介绍的Rosenbrock半隐式Range-Kuta算法，既保持了半隐式Range-Kutta公式的A稳定性，又避免了迭代，是对刚性系统进行实时仿真的较为有效的算法。该算法还有另外一个特性，那就是步长可以根据精度要求进行调节（参[1]p557），从而能在极少的步数下就得到指定时刻系统的状态值的数值解，且满足精度要求。但对于时实仿真来说，固定步长的Rosenbrock算法会更合适，因为时实仿真需要得到各个时刻的精确的系统状态值，若采用变步长，则会遗漏不少需要的观察时刻的系统状态值。
Rosenbrock算法并不复杂，它的描述如下：
[2算法实现]
这里给出的算例如下：
这是一个导弹偏航通道动力学方程描述，有些细节：如系统的参数，变参数特性等，这里不便于公开，敬请谅解。
2．1Matlab实现的关键部分代码：
%--DATA INITIALIZATIONS--
%--调节以下变量以设定不同的状态变量的初始值
anglePSai=0.8; %偏航角
anglePSai1=0.2; %偏航角角速度
angleSigma=0; %弹道偏角
angleDeta=0; %舵偏角
angleBeta=0.8; %侧滑角
aPSaiArr(1)=anglePSai;
aPSai1Arr(1)=anglePSai1;
aSigmaArr(1)=angleSigma;
aDetaArr(1)=angleDeta;
aBetaArr(1)=angleBeta;
y(1)=anglePSai;
y(2)=anglePSai1;
y(3)=angleSigma;
y(4)=angleDeta;
%--MAIN PART OF REAL TIME SIMU --USING ROSENBOCK ALGORITHM , VARIABLE ARGU USE LINEAR INTERPLOTION
i=1;
timeCount=detaT;
arrIndexCount=1;
while timeCount<timeLimit+detaT/2
inb1=b1Arr(arrIndexCount);
inb2=b2Arr(arrIndexCount);
inb3=b3Arr(arrIndexCount);
inc1=c1Arr(arrIndexCount);
inc2=c2Arr(arrIndexCount);
inc3=c3Arr(arrIndexCount);
inF=FArr(arrIndexCount);
inM=MArr(arrIndexCount);
J=[ 0 1 0 0;
-inb2 -inb1 inb2 -inb3;
inc1 0 inc2-inc1 inc3;
C_A0_div_TAO C_A1_div_TAO
fArgus=J;
%---CAL K1---
tmpY=y;
f=(fArgus*tmpY')'+[0 inM -inF 0];
k1=h*inv((eye(matrixSize,matrixSize)-h*b1*J))*f';
%---CAL K2---
tmpY=y+(beta21*k1)';
f=(fArgus*tmpY')'+[0 inM -inF 0];
k2=h*inv((eye(matrixSize,matrixSize)-h*b2*J))*f';
%---CAL Y Vector---
y=y+w1*k1'+w2*k2';
anglePSai=y(1);
anglePSai1=y(2);
angleSigma=y(3);
angleDeta=y(4);
aPSaiArr(arrIndexCount+1)=anglePSai;
aPSai1Arr(arrIndexCount+1)=anglePSai1;
aSigmaArr(arrIndexCount+1)=angleSigma;
aDetaArr(arrIndexCount+1)=angleDeta;
aBetaArr(arrIndexCount+1)=anglePSai-angleSigma;
%--ITERATION CONTROL
arrIndexCount=arrIndexCount+1;
timeCount=timeCount+detaT;
end
2．2C语言实现的关键部分代码：
void SimuStepFrame()
{
double tmpVal;
//USING ROSENBROCK METHOD TO SIMULATION
tmpVal=(timeCount-t1)/(t2-t1);
C1=(dataArea[gI+1][1]-dataArea[gI][1])*tmpVal+dataArea[gI][1];
C2=(dataArea[gI+1][2]-dataArea[gI][2])*tmpVal+dataArea[gI][2];
C3=(dataArea[gI+1][3]-dataArea[gI][3])*tmpVal+dataArea[gI][3];
B1=(dataArea[gI+1][4]-dataArea[gI][4])*tmpVal+dataArea[gI][4];
B2=(dataArea[gI+1][5]-dataArea[gI][5])*tmpVal+dataArea[gI][5];
B3=(dataArea[gI+1][6]-dataArea[gI][6])*tmpVal+dataArea[gI][6];
F=(dataArea[gI+1][7]-dataArea[gI][7])*tmpVal+dataArea[gI][7];
M=(dataArea[gI+1][8]-dataArea[gI][8])*tmpVal+dataArea[gI][8];
Jacobi[0][1]=1;
Jacobi[1][0]=-B2;Jacobi[1][1]=-B1;Jacobi[1][2]=B2 ;Jacobi[1][3]=-B3;
Jacobi[2][0]=C1; Jacobi[2][2]=C2-C1;Jacobi[2][3]=C3;
/*The value below will not change.
Jacobi[3][0]=C_A0_div_TAO;Jacobi[3][1]=C_A1_div_TAO;Jacobi[3][3]= C_MINUS_1_div_TAO;
*/
copyMatrix(Jacobi,fArgus);
//---CAL K1---
copy1DArr(y,tmpY);
calF1(); //f=(fArgus*tmpY')'+[0 inM -inF 0];
calK1(b1); //k1=h*inv((eye(matrixSize,matrixSize)-h*b1*J))*f';
//---CAL K2---
calTmpY(); // tmpY=y+(beta21*k1)';
calF2(); //f=(fArgus*tmpY')'+[0 inM -inF 0];
calK2(b2); //k2=h*inv((eye(matrixSize,matrixSize)-h*b2*J))*f';