0
$\begingroup$

I am trying to calculate the exponent of a 3 x 3 matrix using the formula

$\sum_{i=0}^\infty\frac {A^n}{n!}$

I believe that my error may lay in the scalar division with a factorial or the member function performing the actual exponential calls. After a couple of iterations, I see my output values explode for the simple case of having the matrix in question filled with 1 + i elements.

Output This program will calculat

> Blockquotee the exp of a matrix A
Calculating power of matrix
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
Performing scalar division
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)

Performing Summation
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
(1,1)(1,1)(1,1)
1

Calculating power of matrix
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
Performing scalar division
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)
(0,6)(0,6)(0,6)

Performing Summation
(1,7)(1,7)(1,7)
(1,7)(1,7)(1,7)
(1,7)(1,7)(1,7)
2

Calculating power of matrix
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
Performing scalar division
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)
(-18,18)(-18,18)(-18,18)

Performing Summation
(-17,25)(-17,25)(-17,25)
(-17,25)(-17,25)(-17,25)
(-17,25)(-17,25)(-17,25)
3

Calculating power of matrix
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
Performing scalar division
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)
(-108,0)(-108,0)(-108,0)

Performing Summation
(-125,25)(-125,25)(-125,25)
(-125,25)(-125,25)(-125,25)
(-125,25)(-125,25)(-125,25)
4

Calculating power of matrix
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
Performing scalar division
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)
(-324,-324)(-324,-324)(-324,-324)

Performing Summation
(-449,-299)(-449,-299)(-449,-299)

(-449,-299)(-449,-299)(-449,-299) (-449,-299)(-449,-299)(-449,-299) 5

Calculating power of matrix
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
Performing scalar division
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)
(0,-1944)(0,-1944)(0,-1944)

Performing Summation
(-449,-2243)(-449,-2243)(-449,-2243)
(-449,-2243)(-449,-2243)(-449,-2243)
(-449,-2243)(-449,-2243)(-449,-2243)
6

Calculating power of matrix
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)    
Performing scalar division
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)
(5832,-5832)(5832,-5832)(5832,-5832)

Performing Summation
(5383,-8075)(5383,-8075)(5383,-8075)
(5383,-8075)(5383,-8075)(5383,-8075)
(5383,-8075)(5383,-8075)(5383,-8075)
7

Calculating power of matrix
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)    
Performing scalar division
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)
(34992,0)(34992,0)(34992,0)

Performing Summation
(40375,-8075)(40375,-8075)(40375,-8075)
(40375,-8075)(40375,-8075)(40375,-8075)
(40375,-8075)(40375,-8075)(40375,-8075)
8

Calculating power of matrix
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
Performing scalar division
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)
(104976,104976)(104976,104976)(104976,104976)

Performing Summation    
(145351,96901)(145351,96901)(145351,96901)
(145351,96901)(145351,96901)(145351,96901)
(145351,96901)(145351,96901)(145351,96901)
9

#pragma once
#include <iostream>
#include <complex>
#include <cmath>
#include <cassert>

using namespace std;

class ComplexMatrix
{
    private:
    complex<long double>** Arr;
    int mi = 3;
    int mj = 3;
public:
    ComplexMatrix();
    //~ComplexMatrix();
    ComplexMatrix(int i, int j);
    ComplexMatrix(const ComplexMatrix&);
    void Initialise(complex<long double>);
    void DisplayMatrix();
    void DeleteMatrix();
    void EnterComplexMatrix(int, int);

    ComplexMatrix matrixPower(ComplexMatrix&, int);
    ComplexMatrix ScalarDivision_Fac(ComplexMatrix& , complex<long 
    double>);
    ComplexMatrix MatrixAddition(ComplexMatrix&, ComplexMatrix&);

    ComplexMatrix matrixSq(ComplexMatrix&);
    ComplexMatrix Multiply(ComplexMatrix, ComplexMatrix);

    ComplexMatrix matrixe_A(ComplexMatrix&, int);

    ComplexMatrix operator*(const ComplexMatrix&);
    ComplexMatrix operator=(const ComplexMatrix&);

    friend ComplexMatrix operator+(const ComplexMatrix&, const ComplexMatrix&);
}; 

//Constructor to initialise a default 3 x 3 complex matrix with 0s for real and imaginary values
ComplexMatrix::ComplexMatrix()
{
    Arr = new complex<long double> * [mi];

    for (int x = 0; x < mi; x++)
    {
        Arr[x] = new complex<long double>  [mj];
    }

    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            Arr[x1][y1] = (0.0, 0.0);
        }    
    }     
}

//Constructor to initialise a user defined complex matrix of a particular size with 0s for real                                    
and imaginary values
ComplexMatrix::ComplexMatrix(int i, int j)
{
mi = i;
mj = j;
Arr = new complex<long double> * [i];

for (int x = 0; x < i; x++)
{
    Arr[x] = new complex<long double>[j];
}
for (int x1 = 0; x1 < i; x1++)
{
    for (int y1 = 0; y1 < j; y1++)
    {
        Arr[x1][y1] = (0.0, 0.0);
    }
    }
}

//Copy Constructor
ComplexMatrix::ComplexMatrix(const ComplexMatrix& CM)
//lets only have ARR mean one thing to make it easier to read, understand, 
and to avoid this->         
clutter. 
{
    Arr = new complex<long double> * [mi];

    for (int x = 0; x < mi; x++)
    {
        Arr[x] = new complex<long double>[mj];
    }

for (int x1 = 0; x1 < mi; x1++)
{
    for (int y1 = 0; y1 < mj; y1++)
    {
        Arr[x1][y1] = CM.Arr[x1][y1];
    }
    }
}
/*    
ComplexMatrix::~ComplexMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        delete[] Arr[x];
    }
    delete[] Arr;
}
*/
//Initialise matrix elements to a particular value
void ComplexMatrix::Initialise(complex<long double> x)
{
    for (int i = 0; i < mi; i++)
    {
        for (int j = 0; j < mj; j++)
        {
            Arr[i][j] = x;
        }
    }
}

//Display the matrix member function
void ComplexMatrix::DisplayMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        for (int y = 0; y < mj; y++)
        {
            cout << Arr[x][y];
        }
        cout << endl; 
    }
}

//Delete the memory allocated by the matrix member function
void ComplexMatrix::DeleteMatrix()
{
    for (int x = 0; x < mi; x++)
    {
        delete[] Arr[x];
    }
    delete[] Arr;
}

//Enter complex matrix elements member function
void ComplexMatrix::EnterComplexMatrix(int i, int j)
{
    double real, img;
complex < long double> temp = (0.0, 0.0);
cout << "Your matrix will have " << i * j << " elements" << endl;

//Prompt for user input and assign values for real and imaginary values
for (int x = 0; x < i; x++)
{
    for (int y = 0; y < j; y++)
    {
        cout << "Enter the details for the real part of element" << "[" << 
 x << "]" << "[" << y 
        << "]" << endl;
        cin >> real;
        cout << "Enter the details for the real part of element" << "[" << 
 x << "]" << "[" << y 
        << "]" << endl;
        cin >> img;
        temp = (real, img);
        Arr[x][y] = temp;
        }
    }
}

ComplexMatrix ComplexMatrix::Multiply(ComplexMatrix x, ComplexMatrix y)
{
    ComplexMatrix z(3, 3);

    for (int x1 = 0; x1 < 3; ++x1)
    {
        for (int y1 = 0; y1 < 3; ++y1)
        {
            for (int z1 = 0; z1 < 3; ++z1)
            {
                Arr[x1][y1] += x.Arr[x1][z1] * y.Arr[z1][y1];
            }
        }
    }    
    return z;
}

ComplexMatrix ComplexMatrix::ScalarDivision_Fac(ComplexMatrix& x, 
complex<long double> n)
{
ComplexMatrix newCompArr(3, 3);
complex <long double> fac = 0.0;
int n1 = static_cast <int>(n.real());
n1 = static_cast <int>(n1);
complex <long double> i1;

for (int i = 1; i < n1; i++)
{
    i1 = i;
    fac = fac * i1;
}

for (int x1 = 0; x1 < mi; x1++)
{
    for (int y1 = 0; y1 < mj; y1++)
    {
        newCompArr.Arr[x1][y1] = x.Arr[x1][y1] / fac;
    }
}

return newCompArr;
}

ComplexMatrix ComplexMatrix::matrixSq(ComplexMatrix& x)
{
    ComplexMatrix result(mi, mj);
    result = x * x;
    return result;
}

ComplexMatrix ComplexMatrix::matrixPower(ComplexMatrix& a, int n)
{
    ComplexMatrix result(mi, mj);
    ComplexMatrix temp(mi, mj);

    temp = a;

    if (n % 2 == 0)
    {
        for (int i = 1; i < n / 2; i++)
        {
            result = temp * a;
            temp = result;
        }
        result = temp;
        result = result.matrixSq(result);
    }
    else
    {
        for (int j = 0; j < (n - 1) ; j++)
        {
            result = temp * a;
            temp = result;
        }
        result = temp;
    }
    return result;
}

ComplexMatrix ComplexMatrix::matrixe_A(ComplexMatrix& A, int n)
{
    ComplexMatrix expA(mi, mj);
    ComplexMatrix sum(mi, mj);
    sum.Initialise({ 0.0, 0.0 });
    ComplexMatrix A_n(mi, mj);
    ComplexMatrix A_n_div_n(mi, mj);
    ComplexMatrix temp(mi, mj);
    ComplexMatrix zero(mi, mj);
    zero.Initialise({ 0.0, 0.0 });
    complex <long double> j;

    for (int i = 1; i < n; i++)
    {
        cout << "Calculating power of matrix" << endl;
        A_n = A.matrixPower(A, i);
        A_n.DisplayMatrix();

        A_n_div_n = A_n;
        cout << "Performing scalar division" << endl;
        A_n_div_n.ScalarDivision(A_n_div_n, i);
        A_n_div_n.DisplayMatrix();
        cout << endl;

        temp = zero;
        temp = A_n_div_n;
        cout << "Performing Summation" << endl;
        sum = sum + temp;
        sum.DisplayMatrix();

        cout << i << endl;
        cout << endl;

        temp = zero;
        A_n = A.matrixPower(A, i);
    }
    return sum;
}

ComplexMatrix ComplexMatrix::operator*(const ComplexMatrix& CompArr)
{
    ComplexMatrix newCompArr(3, 3);

    for (int x1 = 0; x1 < 3; ++x1)
    {
    for (int y1 = 0; y1 < 3; ++y1)
        {
            newCompArr.Arr[x1][y1] = {0.0, 0.0};
            for (int z1 = 0; z1 < 3; ++z1)
            {
                newCompArr.Arr[x1][y1] += Arr[x1][z1] * CompArr.Arr[z1][y1];
            }
        }
    }
    return newCompArr;
}

ComplexMatrix ComplexMatrix::operator=(const ComplexMatrix& CM)
{
    mi = 3;
    mj = 3;

    //ComplexMatrix Arr(3, 3);

    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            Arr[x1][y1] = CM.Arr[x1][y1];
        }
    }
    //return Arr;
    return *this;
}

ComplexMatrix operator+(const ComplexMatrix &x, const ComplexMatrix &y)
{
    int mi = 3;
    int mj = 3;
    ComplexMatrix newCompArr(3, 3);
    for (int x1 = 0; x1 < mi; x1++)
    {
        for (int y1 = 0; y1 < mj; y1++)
        {
            newCompArr.Arr[x1][y1] = {  (x.Arr[x1][y1].real() + y.Arr[x1][y1].real()) ,
                                    (x.Arr[x1][y1].imag() + y.Arr[x1][y1].imag()) };
        }
    }
    return newCompArr;
}

#include <iostream>
#include "Header.h"
#include <complex>
#include <cmath>

using namespace std;

int main()
{
    std::cout << "This program will calculate the exp of a matrix A\n";
   complex<long double> x = {1, 1};
    complex<long double> y = {2, 2};
    complex<long double> z = { 0.0, 0.0 };

    ComplexMatrix z1(3, 3);
    ComplexMatrix z2(3, 3);
    ComplexMatrix z3(3, 3);

    z2.Initialise(x);
    z3 = z2.matrixe_A(z2, 10);
    //z3.DisplayMatrix();

    z1.DeleteMatrix();
    z2.DeleteMatrix();
}
Improve this question
$\endgroup$
5
  • $\begingroup$ We are not too keen on debugging code here. However, keep in mind that when computing exponentials with the Taylor series it is very well possible that the intermediate values are much larger than the final result. Do you have a way to compare your intermediate results with a pencil-and-paper computation? $\endgroup$
    – Federico Poloni
    Commented May 26, 2020 at 11:40
  • $\begingroup$ Thanks! Don't know why I didn't do that before with the outputs. Should be able to use some online calculators. $\endgroup$
    – Kishan Bhatt
    Commented May 26, 2020 at 11:54
  • 1
    $\begingroup$ Using the Taylor series of the exponential is likely one of the least suited ways to compute the matrix exponential. There is a well known paper called something like "Ten ways to compute the matrix exponential" that shows how to do what you want to do. $\endgroup$
    – Wolfgang Bangerth
    Commented May 26, 2020 at 15:34
  • $\begingroup$ Actually, 19 ways: cs.cornell.edu/cv/ResearchPDF/19ways+.pdf $\endgroup$
    – Wolfgang Bangerth
    Commented May 26, 2020 at 15:35
  • $\begingroup$ I was just following a textbook exercise. Ill also refer to the paper and maybe implement one of the others as a project. Thanks for the feedback. $\endgroup$
    – Kishan Bhatt
    Commented May 26, 2020 at 15:37

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.