C++

Forum

 
Error in Complex Matrix exponentiation calculation within a class

Shishykish (140)
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)

Header
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
 
 #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 Exp_by_sq(ComplexMatrix&, int);
	void ScalarDivision(complex<long double>);
	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<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;
}

void ComplexMatrix::ScalarDivision(complex<long double> 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
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 = {2, 2};
    complex<long double> 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
Ganado (6783)
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
jonnin (11329)
Taking the exponential of a matrix is not the same thing as repeatedly multiplying the same thing.

thankfully, its actually much faster than this :)
Shishykish (140)
@Ganado

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.
Ganado (6783)
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
jonnin (11329)
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
lastchance (6980)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
 
#include <iostream>
#include <iomanip>
#include <vector>
#include <complex>
#include <cassert>
#include <cmath>
using namespace std;

using complx = complex<double>;
using matrix = vector<vector<complx>>;


//==========================================================


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<complx>( 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<complx>( 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
Shishykish (140)
@Ganado

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.
Shishykish (140)
@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.
Ganado (6783)
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.
Shishykish (140)
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.

Header
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
 
#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 Exp_by_sq(ComplexMatrix&, int);
	void ScalarDivision(complex<long double>);
	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<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;
}

void ComplexMatrix::ScalarDivision(complex<long double> 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
 
#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 = {2, 2};
    complex<long double> 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();
}
Shishykish (140)
This is the best complex matrix app for just squaring matrices
jonnin (11329)
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.
lastchance (6980)
matlab is bloated, clunky and expensive. Just use python if you want to check with library routines:

1
2
3
4
5
6
7
8
9
10
 
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
Shishykish (140)
@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!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 
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;
}
Shishykish (140)
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.

Header
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
 
#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(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(ComplexMatrix& x, complex<long double> 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 <long double> j = 0;

	for (int i = 1; i < n; i++)
	{
		A_n = A.matrixPower(A, i);
		complex <long double> 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;
}


Main
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
#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};

    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();
}
Topic archived. No new replies allowed.