30 # define _USE_MATH_DEFINES
65 mMaxRealOfComplex(0.0),
66 mImagOfMaxComplex(0.0),
67 mFreqOfMaxComplex(0.0),
68 mOscillationIndicator(0.0),
69 mOscillationIndicator_EV(0.0),
70 mBifurcationIndicator_Hopf(0.0),
71 mBifurcationIndicator_Fold(0.0),
72 mBifurcationIndicator_Hopf_BDT(0.0),
73 mBifurcationIndicator_Fold_BDT(0.0),
99 mMaxrealpart(src.mMaxrealpart),
100 mMaximagpart(src.mMaximagpart),
101 mNposreal(src.mNposreal),
102 mNnegreal(src.mNnegreal),
105 mNcplxconj(src.mNcplxconj),
107 mStiffness(src.mStiffness),
108 mHierarchy(src.mHierarchy),
110 mMaxRealOfComplex(src.mMaxRealOfComplex),
111 mImagOfMaxComplex(src.mImagOfMaxComplex),
112 mFreqOfMaxComplex(src.mFreqOfMaxComplex),
113 mOscillationIndicator(src.mOscillationIndicator),
114 mOscillationIndicator_EV(src.mOscillationIndicator_EV),
115 mBifurcationIndicator_Hopf(src.mBifurcationIndicator_Hopf),
116 mBifurcationIndicator_Fold(src.mBifurcationIndicator_Fold),
117 mBifurcationIndicator_Hopf_BDT(src.mBifurcationIndicator_Hopf_BDT),
118 mBifurcationIndicator_Fold_BDT(src.mBifurcationIndicator_Fold_BDT),
120 mResolution(src.mResolution),
257 for (; pA != pAEnd; ++pA, ++pMatrix)
261 if (!finite(*pA) && !isnan(*pA))
432 else if (
mInfo <= mN)
437 else if (
mInfo == mN + 1)
442 else if (
mInfo == mN + 2)
469 size_t *pTo = Pivot.
array();
470 size_t *pFrom = pTo +
mN - 1;
472 for (; pTo < pFrom; ++pTo, --pFrom)
488 for (i = 0; (
C_INT) i < mN; i++)
505 if (fabs(
mR[i]) > resolution)
508 if (
mR[i] >= resolution)
512 if (
mR[i] <= -resolution)
515 if (fabs(
mI[i]) > resolution)
531 if (fabs(
mI[i]) > resolution)
552 if (
mR[0] > fabs(
mR[mN - 1]))
572 maxt = tott = fabs(1 /
mR[mn]);
575 for (i = 1; (
C_INT) i < mN; i++)
578 distt += maxt - fabs(1 /
mR[i]);
579 tott += fabs(1 /
mR[i]);
589 if (fabs(
mI[0]) > resolution)
599 std::complex<C_FLOAT64> tmpcpl = 1.0;
600 unsigned C_INT32 index_min = 0;
603 for (i = 0; (
C_INT) i < mN; i++)
605 tmpcpl *= std::complex<C_FLOAT64>(
mR[i],
mI[i]);
607 if (fabs(mR[i]) < tmpmin)
609 tmpmin = fabs(mR[i]);
618 for (i = 0; (
C_INT) i < mN - 1; i++)
620 tmpcpl *= (std::complex<C_FLOAT64>(
mR[i],
mI[i]) + std::complex<C_FLOAT64>(mR[i + 1],
mI[i + 1]));
630 for (i = 0; (
C_INT) i < mN; i++)
633 tmpcpl *= std::complex<C_FLOAT64>(
mR[i],
mI[i]);
640 for (i = 0; (
C_INT) i < mN; i++)
642 tmp_product *=
mR[i] / (1 - 0.99 * exp(-
mI[i]));
696 os <<
"KINETIC STABILITY ANALYSIS";
698 os <<
"The linear stability analysis based on the eigenvalues" << std::endl;
699 os <<
"of the Jacobian matrix is only valid for steady states." << std::endl;
701 os <<
"Summary:" << std::endl;
709 os <<
"is asymptotically stable";
711 os <<
"'s stability is undetermined";
715 os <<
"," << std::endl;
716 os <<
"transient states in its vicinity have oscillatory components";
719 os <<
"." << std::endl;
722 os <<
"Eigenvalue statistics:" << std::endl;
724 os <<
" Largest real part: ";
725 os << std::setprecision(6) << A.
mMaxrealpart << std::endl;
727 os <<
" Largest absolute imaginary part: ";
728 os << std::setprecision(6) << A.
mMaximagpart << std::endl;
731 os <<
" The complex eigenvalues with the largest real part are: "
735 os.unsetf(std::ios_base::scientific);
736 os.unsetf(std::ios_base::showpoint);
738 os <<
" are purely real" << std::endl;
741 os <<
" are purely imaginary" << std::endl;
744 os <<
" are complex" << std::endl;
747 os <<
" are equal to zero" << std::endl;
750 os <<
" have positive real part" << std::endl;
753 os <<
" have negative real part" << std::endl;
756 os.setf(std::ios_base::showpoint);
758 os <<
" stiffness = " << A.
mStiffness << std::endl;
759 os <<
" time hierarchy = " << A.
mHierarchy << std::endl;
C_FLOAT64 mBifurcationIndicator_Fold
void sortWithPivot(RandomAccessIterator first, RandomAccessIterator last, CVector< size_t > &pivot)
C_FLOAT64 mMaxRealOfComplex
void calcEigenValues(const CMatrix< C_FLOAT64 > &matrix)
C_FLOAT64 mBifurcationIndicator_Hopf_BDT
C_FLOAT64 mFreqOfMaxComplex
virtual size_t numRows() const
C_FLOAT64 mOscillationIndicator_EV
CVector< C_FLOAT64 > mWork
void resize(size_t size, const bool ©=false)
C_FLOAT64 mBifurcationIndicator_Hopf
const CVector< C_FLOAT64 > & getI() const
C_FLOAT64 mOscillationIndicator
C_FLOAT64 mBifurcationIndicator_Fold_BDT
const size_t & getNimag() const
C_FLOAT64 mImagOfMaxComplex
void stabilityAnalysis(const C_FLOAT64 &resolution)
const size_t & getNposreal() const
const CVector< C_FLOAT64 > & getR() const
virtual void resize(size_t rows, size_t cols, const bool ©=false)
const C_FLOAT64 & getHierarchy() const
virtual void print(std::ostream *ostream) const
const C_FLOAT64 & getMaxrealpart() const
const size_t & getNcplxconj() const
const C_FLOAT64 & getMaximagpart() const
const C_FLOAT64 & getStiffness() const
virtual size_t size() const
int dgees_(char *jobvs, char *sort, L_fp select, integer *n, doublereal *a, integer *lda, integer *sdim, doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, integer *lwork, logical *bwork, integer *info)
std::ostream & operator<<(std::ostream &os, const CEigen &A)
CEigen(const std::string &name="NoName", const CCopasiContainer *pParent=NULL)
CCopasiObject * addVectorReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
const size_t & getNnegreal() const
virtual size_t numCols() const
const size_t & getNzero() const
CCopasiObject * addObjectReference(const std::string &name, CType &reference, const unsigned C_INT32 &flag=0)
const size_t & getNreal() const
#define CONSTRUCTOR_TRACE
bool applyPivot(const CVectorCore< size_t > &pivot)