expm.h

Go to the documentation of this file.

/* $Id: expm.h,v 1.8.2.8 2019/05/28 11:13:02 garloff Exp $ */

#ifndef TBCI_SOLVER_EXPM_H
#define TBCI_SOLVER_EXPM_H

#include "tbci/basics.h"
#include "tbci/matrix.h"

NAMESPACE_TBCI

INST(template <  Matrix<T> > class NN friend double norm_1 (const   Matrix<T>&);)
INST(template <BdMatrix<T> > class NN friend double norm_1 (const BdMatrix<T>&);)

template <typename Matrix>
double norm_1(const Matrix& mat)
{
        double max = 0;
        double sum;

        for (unsigned int i=0; i<mat.rows(); ++i) {
                sum = MATH__ fabs(mat.get(i, 0));
                for (unsigned int j=1; j<mat.columns(); ++j) 
                        sum += MATH__ fabs(mat.get(i, j));
                max = MAX(max, sum);
        }
        return max;
}

INST(template <typename T> class TMatrix friend double norm_1 (const TMatrix<T>&);)
template <typename T>
double norm_1 (const TMatrix<T>& mat)
{
        double max = 0;
        double sum;

        for (unsigned int i=0; i<mat.rows(); ++i) {
                sum = MATH__ fabs(mat.get(i, 0));
                for (unsigned int j=1; j<mat.columns(); ++j) 
                        sum += MATH__ fabs(mat.get(i, j));
                max = MAX(max, sum);
        }
        mat.mark_destroy();
        return max;
}

INST(template <   Matrix<T> > class NN friend void expm (const Matrix<T>&, Matrix<T>&);)
//INST(template < BdMatrix<T> > class NN friend void expm (const BdMatrix<T>&, BdMatrix<T>&);)

//******************************************************************
// Exponential Entwicklung nach Pade (siehe MatLab expm1.m)
//******************************************************************
template<typename Matrix>
void expm (const Matrix& exponent, Matrix& E)
{
        // Matrizen kopieren
        Matrix temp(exponent);

        unsigned int dim = exponent.rows();
        int e,p,q,k;

        //------------------------------------------------
        // Einheitsmatrix generieren
        //------------------------------------------------
        Matrix eye(0.0, dim,dim);
        for(unsigned int j=0;j<dim;j++) eye(j,j)=1.0;

        //Die Matrix normieren, so dass norm_1 <1/2
        double z_norm = norm_1(exponent);
        /*double f =*/ frexp(z_norm,&e);
        int    s = MAX(0, e+1);
        temp /= (pow(2.0, s));

        // Naeherung fuer exp(exponent) (Pade)
        Matrix  X(temp);
        Matrix cX(dim,dim);
        double c=0.5;
        Matrix D (eye-temp*c);
        E = eye+temp*c; 
        q = 6;
        p = 1;

        for(k=2;k<=q;k++)
        {
                c     = c*(q-k+1)/(k*(2*q-k+1));
                X     = X*temp;
                cX    = X*c;
                E   += cX;
                if(p)
                   D += cX;
                else
                   D -= cX;
                p = !p;
        }

        //Gleichungssystem loesen
        E = TBCI__ lu_solve(D,E);

  // Entnormieren
        for(k=1;k<=s;k++) E = E*E;
}


INST2(template < Matrix<T>, Matrix<T> > class NN friend \
      void expm2 (const Matrix<T>& exponent, Matrix<T>& erg, double);)
INST2(template < BdMatrix<T>, Matrix<T> > class NN friend \
      void expm2 (const BdMatrix<T>& exponent, Matrix<T>& erg, double);)

//******************************************************************
// Exponential Entwicklung mit Taylorreihe
//******************************************************************
template<typename MatrixIn, typename MatrixOut>
void expm2 (const MatrixIn& exponent, MatrixOut& erg, double ERROR=1e-3)
{
        unsigned int dim=exponent.rows();

        //Die Matrix normieren, so dass norm_1 <1/2
        int e;
        double z_norm=norm_1(exponent);
        /*double f =*/ frexp(z_norm,&e);
        int    s = MAX(0, e+1);
        MatrixIn scaled(exponent);
        double scale_factor = (pow(2.0, s));
        scaled/=scale_factor;

//      STD__ cout << "\nscale_factor = " << scale_factor << STD__ endl;

        //***************************************************************
        // Naeherung fuer exp(exponent) (Taylor-Reihe)
        //***************************************************************
        MatrixOut temp(scaled);
        double error=1;

        // Reihenglieder 0 und 1 (umgekehrte Reihenfolge!)
        erg = temp;
        for (unsigned int j=0; j<dim; j++) erg(j,j) += 1;

        // weitere Reihenglieder        
        int i;
        for(i=2; i<100 && error>ERROR; i++)
        {
          temp = (temp/i)*scaled;
          erg += temp;
          
          error /= i;
          //STD__ cout << "Genauigkeit=" << error << " " << 1/error << STD__ endl;
        }
//      cout<<"Iterationen "<<i<<endl;

        // Entnormieren
        for(int k=1;k<=s;k++) erg = erg*erg;
}



INST2(template <Matrix<T>, Matrix<T> > class NN friend \
      void expm3 (const Matrix<T>& exponent, Matrix<T>& erg, double);)
INST2(template <BdMatrix<T>, Matrix<T> > class NN friend \
      void expm3 (const BdMatrix<T>& exponent, Matrix<T>& erg, double);)
//******************************************************************
// Exponential Entwicklung als Reihennaeherung
//******************************************************************
template<typename MatrixIn, typename MatrixOut>
void expm3 (const MatrixIn& exponent, MatrixOut& erg, double ERROR=1e-3)
{
        unsigned int dim=exponent.rows();
        double z_norm=0;

        //Die Matrix normieren auf 0.1
        for (unsigned int i=0; i<dim; i++)
                for (unsigned int j=0; j<dim; j++) 
                        z_norm=MAX(z_norm, (double)MATH__ fabs(exponent(i,j)));

        z_norm /= 0.1;
  
        int e;
        /*double f =*/ frexp(z_norm,&e);
        int    s = MAX(0, e+1);
        double scale_factor = (pow(2.0, s));
        double error = 10;
 
        MatrixOut temp(exponent);

//      STD__ cout << z_norm << " " << s << " " << scale_factor << STD__ endl;
        
        for(int iter=0; iter<100 && error>ERROR; iter+=4)
        {

//        STD__ cout << "s=" << s << STD__ endl;
          scale_factor = (pow(2.0, s));
          erg = exponent/(scale_factor);
          for(unsigned int j=0; j<dim; j++) erg(j,j) += 1;
          for(int k=1;k<=s;k++) 
                erg = erg*erg;
          error = norm_1 (erg - temp);
    
          STD__ cout << error << STD__ endl;
          temp=erg;
          s++;
        }

}

NAMESPACE_END

#endif  /* TBCI_SOLVER_EXPM_H */