Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "SIBS.H"
00027
00028
00029
00030 void Foam::SIBS::SIMPR
00031 (
00032 const ODE& ode,
00033 const scalar xStart,
00034 const scalarField& y,
00035 const scalarField& dydx,
00036 const scalarField& dfdx,
00037 const scalarSquareMatrix& dfdy,
00038 const scalar deltaX,
00039 const label nSteps,
00040 scalarField& yEnd
00041 ) const
00042 {
00043 scalar h = deltaX/nSteps;
00044
00045 scalarSquareMatrix a(n_);
00046 for (register label i=0; i<n_; i++)
00047 {
00048 for (register label j=0; j<n_; j++)
00049 {
00050 a[i][j] = -h*dfdy[i][j];
00051 }
00052 ++a[i][i];
00053 }
00054
00055 labelList pivotIndices(n_);
00056 LUDecompose(a, pivotIndices);
00057
00058 for (register label i=0; i<n_; i++)
00059 {
00060 yEnd[i] = h*(dydx[i] + h*dfdx[i]);
00061 }
00062
00063 LUBacksubstitute(a, pivotIndices, yEnd);
00064
00065 scalarField del(yEnd);
00066 scalarField ytemp(n_);
00067
00068 for (register label i=0; i<n_; i++)
00069 {
00070 ytemp[i] = y[i] + del[i];
00071 }
00072
00073 scalar x = xStart + h;
00074
00075 ode.derivatives(x, ytemp, yEnd);
00076
00077 for (register label nn=2; nn<=nSteps; nn++)
00078 {
00079 for (register label i=0; i<n_; i++)
00080 {
00081 yEnd[i] = h*yEnd[i] - del[i];
00082 }
00083
00084 LUBacksubstitute(a, pivotIndices, yEnd);
00085
00086 for (register label i=0; i<n_; i++)
00087 {
00088 ytemp[i] += (del[i] += 2.0*yEnd[i]);
00089 }
00090
00091 x += h;
00092
00093 ode.derivatives(x, ytemp, yEnd);
00094 }
00095 for (register label i=0; i<n_; i++)
00096 {
00097 yEnd[i] = h*yEnd[i] - del[i];
00098 }
00099
00100 LUBacksubstitute(a, pivotIndices, yEnd);
00101
00102 for (register label i=0; i<n_; i++)
00103 {
00104 yEnd[i] += ytemp[i];
00105 }
00106 }
00107
00108
00109