### Error in Complex Matrix exponentiation calculation within a class

Hello

I am writing a program to calculate the exp of a matrix for a given formula. I can successfully calculate the scalar division and multiplication of the matrix as an objects within a class called ComplexMatrix with complex elements. I understand the algorithm for calculating the matrix raised by a power.

I used a test matrix of 3 x 3 with 2+2i and expect to get 24i as a result from some online calculators. When I call the function I get.

(2,2)(2,2)(2,2)
(2,2)(2,2)(2,2)
(2,2)(2,2)(2,2)

(2,26)(2,26)(2,26)
(2,26)(2,26)(2,26)
(2,26)(2,26)(2,26)

 ``123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261`` `````` #pragma once #include #include #include #include using namespace std; class ComplexMatrix { private: complex** Arr; int mi = 3; int mj = 3; public: ComplexMatrix(); //~ComplexMatrix(); ComplexMatrix(int i, int j); ComplexMatrix(const ComplexMatrix&); void Initialise(complex); void DisplayMatrix(); void DeleteMatrix(); void EnterComplexMatrix(int, int); ComplexMatrix Exp_by_sq(ComplexMatrix&, int); void ScalarDivision(complex); ComplexMatrix Multiply(ComplexMatrix, ComplexMatrix); ComplexMatrix operator*(const ComplexMatrix&); ComplexMatrix operator=(const ComplexMatrix&); }; //Constructor to initialise a default 3 x 3 complex matrix with 0s for real and imaginary values ComplexMatrix::ComplexMatrix() { Arr = new complex * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex [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 * [i]; for (int x = 0; x < i; x++) { Arr[x] = new complex[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 * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex[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 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; } void ComplexMatrix::ScalarDivision(complex n) { ComplexMatrix newCompArr(3, 3); for (int x1 = 0; x1 < mi; x1++) { for (int y1 = 0; y1 < mj; y1++) { Arr[x1][y1] = Arr[x1][y1] / n; } } } ComplexMatrix ComplexMatrix::Exp_by_sq(ComplexMatrix& x, int n) { ComplexMatrix y(mi, mj); if (n == 1) { for (int i = 0; i < mi; i++) { for (int j = 0; j < mj; j++) { y.Arr[i][j] = (1.0, 1.0); } } } else if (n % 2 == 0) { for (int x1 = 0; x1 < (n / 2); x1++) { y = x.Multiply(x, x); } } else { for (int y1 = 0; y1 < ((n - 1) / 2); y1++) { y = x.Multiply(x, x); } y = x.Multiply(x, x); } return y; } ComplexMatrix ComplexMatrix::operator*(const ComplexMatrix& CompArr) { ComplexMatrix newCompArr(3, 3); for (int x1 = 0; x1 < 3; ++x1) { for (int y1 = 0; y1 < 3; ++y1) { 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& CompArr) { mi = 3; mj = 3; /* assert(mi != mj); { cout << "There are more than 3 rows" << endl; } assert(mj != mi); { cout << "There are more than 3 rows" << endl; } */ ComplexMatrix newCompArr(3, 3); for (int x1 = 0; x1 < mi; x1++) { for (int y1 = 0; y1 < mj; y1++) { newCompArr.Arr[x1][y1] = CompArr.Arr[x1][y1]; } } return *this; } ``````

Main
 ``1234567891011121314151617181920`` ``````include "Header.h" #include #include using namespace std; int main() { std::cout << "This program will calculate the exp of a matrix A\n"; complex x = {2, 2}; complex y = {1, 1}; ComplexMatrix z1(3, 3); z1.Initialise(x); z1.DisplayMatrix(); z1.Exp_by_sq(z1, 2); z1.DisplayMatrix(); z1.DeleteMatrix(); } ``````
Last edited on
`y.Arr[i][j] = (1.0, 1.0);`
Misuse of the comma operator. This is not making a pair or a complex number, it's just assigning it the value 1.0.
Turn on compiler warnings,
 ` warning: left operand of comma operator has no effect [-Wunused-value]`

But... okay, you're trying to calculate the matrix exponential. I don't exactly understand where you're calculating it. Taking the exponential of a matrix is not the same thing as repeatedly multiplying the same thing.
https://en.wikipedia.org/wiki/Matrix_exponential#Computing_the_matrix_exponential

Honestly I just let MATLAB's expm function do all the work for me :^)
MATLAB (and Octave, etc.) apparently use the PadÃ© approximation with scaling and squaring.
Last edited on
Taking the exponential of a matrix is not the same thing as repeatedly multiplying the same thing.

thankfully, its actually much faster than this :)

From that article, the formula I am using is the same as in the projection heading. Just without the factorial. I wanted a class member function that would raise the matrix to a particular power. This using that in conjuction with my scalar division function and another function to perform a summation.

Matlab? I know it does these things. I am trying to teach myself scientific computing and C++ was supposedley good for this. Is matlab useful too?

 Taking the exponential of a matrix is not the same thing as repeatedly multiplying the same thing. thankfully, its actually much faster than this :)

Ha! I know that this I am repeatedly multiplying here. I am just working to a rough formula. I needed to raise the matrix to a power, I was using the methodology of https://www.rookieslab.com/posts/fast-power-algorithm-exponentiation-by-squaring-cpp-python-implementation to do so.
I think we're talking about two different things...?
Repeated multiplication (A * A * A *A * ...) is not the same thing as taking the exponential of a matrix (e^A), although you might use repeated exponentiation as part of finding the exponential of a matrix.
Last edited on
Matlab? I know it does these things. I am trying to teach myself scientific computing and C++ was supposedley good for this. Is matlab useful too?

matlab is a very expensive matrix engine with a lot more. It has its own simple but powerful programming language and it is (or was) very slow; the library I wrote actually was a translator for .m files into c++ to take what the math guys did and get it done fast.

Matlab is good to use what has been done for you to do new things. Its not awesome at doing large things quickly, so we used it to prove out the math before coding it in c++. It also let us run the same problem in both to validate the c++. This is a big deal, debugging aside matlab does heavy numerical correctness inside and that is slower but more accurate results. If we had a numerical problem in our code, we could spot it this way. My code was lax on that in favor of speed, and a few times I had to go back to add that in for specific routines.
Last edited on
 ``123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153`` ``````#include #include #include #include #include #include using namespace std; using complx = complex; using matrix = vector>; //========================================================== void print( const matrix &M ) { for ( auto &row : M ) { for ( auto z : row ) cout << setw( 20 ) << z << " "; cout << '\n'; } } //========================================================== matrix multiply( const matrix& A, const matrix &B ) { assert( A[0].size() == B.size() ); matrix M( A.size(), vector( B[0].size(), 0.0 ) ); for ( int i = 0; i < A.size(); i++ ) { for ( int j = 0; j < B[0].size(); j++ ) { M[i][j] = 0.0; for ( int k = 0; k < A[0].size(); k++ ) M[i][j] += A[i][k] * B[k][j]; } } return M; } //========================================================== matrix add( const matrix& A, const matrix &B ) { assert( A .size() == B .size() ); assert( A[0].size() == B[0].size() ); auto M = A; for ( int i = 0; i < A.size(); i++ ) { for ( int j = 0; j < A[0].size(); j++ ) M[i][j] += B[i][j]; } return M; } //========================================================== matrix multScalar( complx z, const matrix& A ) { auto M = A; for ( int i = 0; i < A.size(); i++ ) { for ( int j = 0; j < A[0].size(); j++ ) M[i][j] *= z; } return M; } //========================================================== double L2norm( const matrix& A ) { double norm = 0.0; for ( int i = 0; i < A.size(); i++ ) { for ( int j = 0; j < A[0].size(); j++ ) norm += abs( A[i][j] * A[i][j] ); } return sqrt( norm ); } //========================================================== matrix identity( int n ) { matrix I( n, vector( n, 0.0 ) ); for ( int i = 0; i < n; i++ ) I[i][i] = complx( 1.0 ); return I; } //========================================================== matrix exp( const matrix& A ) { const double eps = 1.0e-10; assert( A.size() == A[0].size() ); matrix term = identity( A.size() ); matrix sum = term; for ( int n = 1; L2norm( term ) > eps; n++ ) { term = multiply( A, term ); term = multScalar( complx( 1.0 / n ), term ); sum = add( sum, term ); } return sum; } //========================================================== matrix pow( matrix A, int n ) { matrix result = identity( A.size() ); while ( n ) { if ( n & 1 ) result = multiply( result, A ); A = multiply( A, A ); n >>= 1; } return result; } //========================================================== int main() { complx z( 2.0, 2.0 ); matrix M = { { z, z, z }, { z, z, z }, { z, z, z } }; cout << "M:\n" ; print( M ); cout << "\nexp(M):\n"; print( exp( M ) ); cout << "\nM^2:\n" ; print( pow( M, 2 ) ); cout << "\nM^10:\n" ; print( pow( M, 10 ) ); }``````

 ```M: (2,2) (2,2) (2,2) (2,2) (2,2) (2,2) (2,2) (2,2) (2,2) exp(M): (129.787,-37.5748) (128.787,-37.5748) (128.787,-37.5748) (128.787,-37.5748) (129.787,-37.5748) (128.787,-37.5748) (128.787,-37.5748) (128.787,-37.5748) (129.787,-37.5748) M^2: (0,24) (0,24) (0,24) (0,24) (0,24) (0,24) (0,24) (0,24) (0,24) M^10: (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) (0,6.44973e+08) ```
Last edited on

 Repeated multiplication (A * A * A *A * ...) is not the same thing as taking the exponential of a matrix (e^A), although you might use repeated exponentiation as part of finding the exponential of a matrix.

I was using the repeated exponentiation to do exactly that! the equation has this in.

@Jonnin
 Matlab is good to use what has been done for you to do new things. Its not awesome at doing large things quickly, so we used it to prove out the math before coding it in c++. It also let us run the same problem in both to validate the c++. This is a big deal, debugging aside matlab does heavy numerical correctness inside and that is slower but more accurate results. If we had a numerical problem in our code, we could spot it this way. My code was lax on that in favor of speed, and a few times I had to go back to add that in for specific routines

Is matlab written in a similar language? Sounds like a great to test models etc. I might check that out as well.
@Last Chance

This looks great! Thanks very much. Some of the vector stuff is unfamiliar to me, but I can look those up. I love it that you didn't have to code "delete[] Matrix;" every other line! I will try developing my own class first and this will give me something to work towards.
MATLAB is its own IDE application and language. It can also interface with C/C++ for speed purposes. But yeah, I feel the same as jonnin. It can be quicker to type something small up in MATLAB for purposes of prototyping (e.g. has matrix exponential built in), but if you're using C++ you might as well continue using C++, but you can verify things with MATLAB. There's also free alternatives, like Octave.
Hello,

So my current code looks like this. In line 21 of main, when I call the exp_by_sq member function on z2 raised to power 2, I get the correct answer of 24i. Here is what I should expect for other powers for a complex matrix 3 x 3 raised to power 3,4 and 5

power = 3 expected = -144 + 144i actual = 0 + 48i
power = 4 expected = -1728 actual = 0 + 48i
power = 5 expected = 10 369 - 10 368i actual = 0 + 72i

Not really sure where I am going wrong with the if-else loop here.

 ``123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282`` ``````#pragma once #include #include #include #include using namespace std; class ComplexMatrix { private: complex** Arr; int mi = 3; int mj = 3; public: ComplexMatrix(); //~ComplexMatrix(); ComplexMatrix(int i, int j); ComplexMatrix(const ComplexMatrix&); void Initialise(complex); void DisplayMatrix(); void DeleteMatrix(); void EnterComplexMatrix(int, int); ComplexMatrix Exp_by_sq(ComplexMatrix&, int); void ScalarDivision(complex); ComplexMatrix MatrixAddition(ComplexMatrix&, ComplexMatrix&); ComplexMatrix Multiply(ComplexMatrix, ComplexMatrix); 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 * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex [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 * [i]; for (int x = 0; x < i; x++) { Arr[x] = new complex[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 * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex[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 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; } void ComplexMatrix::ScalarDivision(complex n) { ComplexMatrix newCompArr(3, 3); for (int x1 = 0; x1 < mi; x1++) { for (int y1 = 0; y1 < mj; y1++) { Arr[x1][y1] = Arr[x1][y1] / n; } } } ComplexMatrix ComplexMatrix::MatrixAddition(ComplexMatrix &x, ComplexMatrix &y) { ComplexMatrix z(3, 3); for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Arr[i][j] =x.Arr[i][j] + y.Arr[i][j]; } } return z; } ComplexMatrix ComplexMatrix::Exp_by_sq(ComplexMatrix& x, int n) { ComplexMatrix result(mi, mj), base(mi, mj); if (n == 1) { for (int i = 0; i < mi; i++) { for (int j = 0; j < mj; j++) { result.Arr[i][j] = (1.0, 1.0); } } } else if ((n % 2) == 0) { for (int x1 = 0; x1 < (n / 2); x1++) { result = result + (x * x); } } else { for (int y1 = 0; y1 < ((n - 1) / 2); y1++) { result = result + (x * x); } result = result + (x * x); } return result; } ComplexMatrix ComplexMatrix::operator*(const ComplexMatrix& CompArr) { ComplexMatrix newCompArr(3, 3); for (int x1 = 0; x1 < 3; ++x1) { for (int y1 = 0; y1 < 3; ++y1) { 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; } ``````

Main
 ``123456789101112131415161718192021222324252627`` ``````#include #include "Header.h" #include #include using namespace std; int main() { std::cout << "This program will calculate the exp of a matrix A\n"; complex x = {2, 2}; complex y = {1, 1}; ComplexMatrix z1(3, 3); ComplexMatrix z2(3, 3); ComplexMatrix z3(3, 3); ComplexMatrix z4(3, 3); z2.Initialise(x); z3 = z2.Exp_by_sq(z2, 2); z3.DisplayMatrix(); z1.DeleteMatrix(); z2.DeleteMatrix(); z3.DeleteMatrix(); z4.DeleteMatrix(); } ``````
This is the best complex matrix app for just squaring matrices
Is matlab written in a similar language? Sounds like a great to test models etc. I might check that out as well

matlab scrips are more like basic than c++.
the variables have no type, so a variable x may be a string like "hello" then 5 lines later it may be a double like 3.14 and 2 lines after that it could be a 3-d matrix of complex<double> and so on. Also a typo, like var1 varl or something will just create a new variable that looks similar but isnt quite right, creating tons of debugging fun.
on the flipside, its powerful: you can draw 2 and 3-d graphs in a line or two, load and save to a file easily, and there are tons of add-ons beyond linear algebra including contols, stastics, and other common math problems that use applied linear algebra or even pure math. My only issue with it was that its slower than c++. depending on what you are doing this varies from a little slower to a whole lot slower. Well, that and the price tag. They are very 'proud' of their product.
matlab is bloated, clunky and expensive. Just use python if you want to check with library routines:

 ``12345678910`` ``````import numpy as np from numpy.linalg import matrix_power from scipy.linalg import expm A = np.array( [ [ 2+2j, 2+2j, 2+2j ], [ 2+2j, 2+2j, 2+2j ], [ 2+2j, 2+2j, 2+2j ] ] ) print( "A:\n", A ) print( "A^2:\n", matrix_power( A, 2 ) ) print( "A^10:\n", matrix_power( A, 10 ) ) print( "exp(A):\n", expm( A ) )``````

 ```A: [[2.+2.j 2.+2.j 2.+2.j] [2.+2.j 2.+2.j 2.+2.j] [2.+2.j 2.+2.j 2.+2.j]] A^2: [[0.+24.j 0.+24.j 0.+24.j] [0.+24.j 0.+24.j 0.+24.j] [0.+24.j 0.+24.j 0.+24.j]] A^10: [[0.+6.44972544e+08j 0.+6.44972544e+08j 0.+6.44972544e+08j] [0.+6.44972544e+08j 0.+6.44972544e+08j 0.+6.44972544e+08j] [0.+6.44972544e+08j 0.+6.44972544e+08j 0.+6.44972544e+08j]] exp(A): [[129.7867801-37.57475244j 128.7867801-37.57475244j 128.7867801-37.57475244j] [128.7867801-37.57475244j 129.7867801-37.57475244j 128.7867801-37.57475244j] [128.7867801-37.57475244j 128.7867801-37.57475244j 129.7867801-37.57475244j]]```
Last edited on
@LastChance

Funnily enough, I am currently learning python alongside C++. I will check out the matrix library with it as well. C++, python and rescuing stray cats seems to be my routine in the time of lockdown. I will double check my results against python instead. I am currently trying get my head around what I am doing wrong here. It's 2349 over here. I think I will sleep on it.

Thanks very much for your help!

 ``1234567891011121314151617181920212223242526272829303132`` ``````ComplexMatrix ComplexMatrix::Exp_by_sq(ComplexMatrix& x, int n) { ComplexMatrix result(mi, mj), base(mi, mj); if (n == 1) { for (int i = 0; i < mi; i++) { for (int j = 0; j < mj; j++) { result.Arr[i][j] = (1.0, 1.0); } } } else if ((n % 2) == 0) { for (int x1 = 0; x1 < (n / 2); x1++) { result = result + (x * x); } } else { for (int y1 = 0; y1 < ((n - 1) / 2); y1++) { result = result + (x * x); } result = result + (x * x); } return result; }``````
Hi,

I thankfully sorted the raising the matrix to a particular power problem. Now just to work on the implementation for the equation below which is an approximation I was working towards

exp^A = (sum from n to inf) (Matrix A^n / n)

Where n would be some suitably large number.

 ``123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315`` ``````#pragma once #include #include #include #include using namespace std; class ComplexMatrix { private: complex** Arr; int mi = 3; int mj = 3; public: ComplexMatrix(); //~ComplexMatrix(); ComplexMatrix(int i, int j); ComplexMatrix(const ComplexMatrix&); void Initialise(complex); void DisplayMatrix(); void DeleteMatrix(); void EnterComplexMatrix(int, int); ComplexMatrix matrixPower(ComplexMatrix&, int); ComplexMatrix ScalarDivision(ComplexMatrix& , complex); 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 * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex [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 * [i]; for (int x = 0; x < i; x++) { Arr[x] = new complex[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 * [mi]; for (int x = 0; x < mi; x++) { Arr[x] = new complex[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 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(ComplexMatrix& x, complex n) { 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] / n; } } return newCompArr; } ComplexMatrix ComplexMatrix::MatrixAddition(ComplexMatrix &x, ComplexMatrix &y) { ComplexMatrix z(3, 3); for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { Arr[i][j] =x.Arr[i][j] + y.Arr[i][j]; } } return z; } 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); complex j = 0; for (int i = 1; i < n; i++) { A_n = A.matrixPower(A, i); complex j = i; A_n_div_n = A_n.ScalarDivision(A_n, j); temp = A_n_div_n; sum = sum + temp; } 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; }``````
 ``1234567891011121314151617181920212223242526272829`` ``````#include #include "Header.h" #include #include using namespace std; int main() { std::cout << "This program will calculate the exp of a matrix A\n"; complex x = {1, 1}; complex y = {2, 2}; ComplexMatrix z1(3, 3); ComplexMatrix z2(3, 3); ComplexMatrix z3(3, 3); ComplexMatrix z4(3, 3); z1.Initialise(x); z2.Initialise(y); z3 = z1.matrixe_A(z1, 100); z3.DisplayMatrix(); z1.DeleteMatrix(); z2.DeleteMatrix(); z3.DeleteMatrix(); z4.DeleteMatrix(); } ``````