C++

Forum

 
matrix power

mrsh (27)
i tried to run matrix power,but i dont get the correct answer for my matrix,this is my code

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
 
#include <fstream.>
#include <iostream>
#include <conio.h>
#include <iomanip>
#include <math.h>
#define n 3

using namespace std;
void main()
{

	double A[n+1][n+1], z[n + 1][n + 1]; 
	int i, j;
	cout.setf(ios::fixed);
	cout.precision(4);

	A[1][1] = 1.0000; A[1][2] = 0.5000;   A[1][3] = 0.1250;
	A[2][1] = 0.0000; A[2][2] = 0.8750;   A[2][3] = 0.2500;
	A[3][1] = 0.0000; A[3][2] = 0.0000;   A[3][3] = 0.1250;


	


	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= n; j++)
		{

			
			cout << "    " << A[i][j];
		}

		cout << endl;
	}
	cout << endl;

	
	
		
		
     for (i = 1; i <=n; i++)
		{
			for (j = 1; j <=n; j++)
			{
				z[i][j] = pow(A[i][j],2);
					
				cout << "    " << z[i][j];			
			}
			
			cout << endl;
		}
	
	_getch();

}
Write your question here.

closed account (48T7M4Gy)
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
 
//#include <fstream>
#include <iostream>
//#include <conio.h>
//#include <iomanip>
#include <cmath>

#define n 3

using namespace std;

int main()
{ 
    double A[n+1][n+1], z[n + 1][n + 1]; 
    
    cout.setf(ios::fixed);
	cout.precision(4);

	A[1][1] = 1.0000; A[1][2] = 0.5000;   A[1][3] = 0.1250;
	A[2][1] = 0.0000; A[2][2] = 0.8750;   A[2][3] = 0.2500;
	A[3][1] = 0.0000; A[3][2] = 0.0000;   A[3][3] = 0.1250;
	
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
		    cout << "    " << A[i][j];
		}
		cout << endl;
	}
	cout << endl;
	
	for (int i = 1; i <=n; i++)
		{
			for (int j = 1; j <=n; j++)
			{
				z[i][j] = pow(A[i][j],2);
					
				cout << "    " << z[i][j];			
			}
			cout << endl;
		}
	//	_getch();

}


There's a bit of a tidy up to get away from devc++ or whatever.

I think the problem is matrix powers are not just each element being raised to that particular power. What you have to do is program multiplying matrices. ie A = B * C by matrix multiplication so, A^2 = A * A, A^3 = A * A * A etc., where each of A, B and C are matrices, not scalars.

If you need to know how to multiply two (conforming matrices) I suggest Google because there are voluminous entries.

BTW Only square matrices can be raised to a power >= 2.

mrsh (27)
thank you,thats mean i have to do a lot of multiplication between matrices?it is possible to do multiplication between 2 matrices like i did below

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
 
#include <fstream.>
#include <iostream>
#include <conio.h>
#include <iomanip>
#include <math.h>
#define n 3

using namespace std;
void main()
{

	double A[n+1][n+1], z[n + 1][n + 1]; 
	int i, j,k;
	cout.setf(ios::fixed);
	cout.precision(4);

	A[1][1] = 1.0000; A[1][2] = 0.5000;   A[1][3] = 0.1250;
	A[2][1] = 0.0000; A[2][2] = 0.8750;   A[2][3] = 0.2500;
	A[3][1] = 0.0000; A[3][2] = 0.0000;   A[3][3] = 0.1250;


	


	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= n; j++)
		{

			
			cout << "    " << A[i][j];
		}

		cout << endl;
	}
	cout << endl;

	
	
		
		
     for (i = 1; i <=n; i++)
		{
			for (j = 1; j <=n; j++)
			{
				z[i][j] = 0;
				for (k = 1; k <= n; k++)
				z[i][j]+= A[i][k]*A[k][j];
				cout << " " << z[i][j];		
			}
			
			cout << endl;
		}
	
	_getch();

}


how about matrix power of 3,4 n so on? do you have any idea?
closed account (48T7M4Gy)
If you adapt what you have and use operator overloading to overload * then you are on the right track.

Some matrices lend themselves to powers but you'll have to surf around or wait for specific expertise on that because I'd have to do the surfing myself. Wolfram might be a start.

I haven't checked your program. You can do that yourself with some test data. I don't think it is correct but I only had a quick look.

Here is an overloaded *operator for matrices.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
//MULTIPLICATION
Matrix operator* ( const Matrix& m1, const Matrix& m2 )
{
    if (m1.m_nCols != m2.m_nRows) { throw std::invalid_argument("incompatible rows & cols"); }
    
    Matrix temp(m1.m_nRows, m2.m_nCols);
    
    for (int i = 0; i <  m1.m_nRows; i++)
    {
        for (int j = 0; j < m2.m_nCols; j++)
        {
            temp.m_dValue[i][j] = 0;
            for (int k = 0; k < m2.m_nRows; k++)
            {
                temp.m_dValue[i][j] += m1.m_dValue[i][k] * m2.m_dValue[k][j];
            }
        }
    }
    return temp;
}
lastchance (6980)
Your matrix multiply is correct (I think).

For higher powers you will need temporary storage for the matrix product to calculate higher powers or you will be over-writing elements of the matrix involved in the multiply.

I suggest that you write a separate void function to multiply two matrices and also one to copy matrices.


To get Z = A^n, where A and Z are matrices; in pseudocode:

- copy A to Z; Z then holds A^1

- loop round a further n-1 times to get successively higher powers:
{
set temp to the matrix product of A and Z
copy temp into Z - This will then hold the latest computed power of A
}


It will definitely be easier if you write separate functions to do the multiplying (and, probably, the copying too.)


If your matrix is upper-triangular (as in your example) then the diagonal elements will individually be the requisite powers. This gives a little bit of a check.
Last edited on
jonnin (11429)
there are several ways to compute this without doing a full multiply for each power, which is terribly slow for bigger matrices.

Also if you transpose before multiply, you can just iterate over both matrices in a much more efficient inner loop (also useful for larger problems).

Not important for 3x3s and school problems, but if you are interested in the area, there are a lot of really cool algorithms out there for linear algebra.
closed account (48T7M4Gy)
@OP You'll be pleased to know I get the same result with my algorithm as you do with yours.

The advantage with operator overloading is when you want to get a power of a (square) matrix A you effectively just write B = A*A* ... *A; No temps, explicit copying or whatever are involved.

It's more or less right that there are algorithms around that will do these jobs faster. With this you end up resorting to eigenvectors and the like but you have to weigh up what sort of powers and what size matrices you plan on handling.

Cholesky and trimatrix decomposition are the common ones instead of Gauss and Gauss-Jordan reduction for inverting matrices and many others taking advantage of matrix symmetries. However when a reasonably fast machine will solve a 500 x 500 augmented matrix in 0.26 second (1000 x 1000 @ 1.9 sec) it's hardly worth the trouble going beyond basics.
mrsh (27)
many thanks for your opinion,i found one in the internet and its work for me, but now i want to divide every matrix by factorials,actually i tried to compute exponential matrices using taylor series and here is my factorials list

1
2
3
4
 
Factorial[2]=2.0000
Factorial[3]=6.0000
Factorial[4]=24.0000
Factorial[5]=120.0000

and here is my matrices power list

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
 1.0000  0.9375  0.2656
  0.0000  0.7656  0.2500
  0.0000  0.0000  0.0156

  1.0000  1.3203  0.3926
  0.0000  0.6699  0.2227
  0.0000  0.0000  0.0020

  1.0000  1.6553  0.5042
  0.0000  0.5862  0.1953
  0.0000  0.0000  0.0002

  1.0000  1.9484  0.6018
  0.0000  0.5129  0.1710
  0.0000  0.0000  0.0000


but what i get after divided by factorials are

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
  0.5000  0.4688  0.1328
  0.0000  0.3828  0.1250
  0.0000  0.0000  0.0078

  0.2500  0.3301  0.0981
  0.0000  0.1675  0.0557
  0.0000  0.0000  0.0005

  0.1250  0.2069  0.0630
  0.0000  0.0733  0.0244
  0.0000  0.0000  0.0000

  0.0625  0.1218  0.0376
  0.0000  0.0321  0.0107
  0.0000  0.0000  0.0000


the result should be

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
0.5000  0.4688  0.1328
  0.0000  0.3828  0.1250
  0.0000  0.0000  0.0078

  0.1667  0.2201  0.0654
  0.0000  0.1117  0.0371
  0.0000  0.0000  0.0003

  0.0417  0.0690  0.0210
  0.0000  0.0244  0.0081
  0.0000  0.0000  0.0000

  0.0083  0.0162  0.0050
  0.0000  0.0043  0.0014
  0.0000  0.0000  0.0000


here is my full code

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
 
#include <fstream.>
#include <iostream>
#include <conio.h>
#include <iomanip>
#include <math.h>
#define n 3
#define p 5

using namespace std;

void main()
{

	double A[n + 1][n + 1], z[n + 1][n + 1], D[n + 1][n + 1],factorial=1;
	int i, j, k, c, d;
	cout.setf(ios::fixed);
	cout.precision(4);

	A[1][1] = 1.0000; A[1][2] = 0.5000;   A[1][3] = 0.1250;
	A[2][1] = 0.0000; A[2][2] = 0.8750;   A[2][3] = 0.2500;
	A[3][1] = 0.0000; A[3][2] = 0.0000;   A[3][3] = 0.1250;

	//=================matrix order 1===================

	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= n; j++)
		{
			D[i][j] = A[i][j];
			cout << "  " << D[i][j];
		}

		cout << endl;
	}
	cout << endl;


	//=================== factorial ====================

	for ( j = 2; j <= p; j++)
	{
		factorial = factorial*j;
		cout << "factorial [" << j << "]" << factorial;
		cout << endl;
	}
	cout << endl;

	//============== matrix order n ================

	

	for (i = 1; i <p; i++)
	{
		
	//==============calculation for matrix multiplication======
	
			for (int c = 1; c <= n; c++)
			{
				for (int d = 1; d <= n; d++)
				{
					z[c][d] = 0;
					for (int k = 1; k <= n; k++)
					{
						//z[c][d] += (D[c][k] * A[k][d])/(factorial*j);
						z[c][d] += D[c][k] * A[k][d];
					}


				}
			}

	//===============to display the matrix===============
		
				for (int c = 1; c <= n; c++)
				{
					for (int d = 1; d <= n; d++)
					{
						D[c][d] = z[c][d];//---------------matrix power
						//D[c][d] = (z[c][d]/factorial*j)/0.1;//----------------matrix power divided by factorial
						cout << "  " << D[c][d];
					
					}
					cout << endl;	
					
				}
				cout << endl;
				
      }
	


	_getch();

}


lastchance (6980)
Yes, but presumably there is one factorial for each power - whereas you appear to be only computing 5! and not all the others.

Taylor series (actually Maclaurin series for exp(A) here):
1 + A + A2/2! + A3/3! + A4/4! + A5/5! + ... = exp(A)
(I am guessing that you might want a sum of this series as well.)

The other problem with factorials is that they grow very fast - either store 1.0/n! or, better, combine them with your successive powers of A:
A2/2! = (A/1!) * A/2
A3/3! = (A2/2!) * A/3
A4/4! = (A3/3!) * A/4
etc.

You can keep An/n! in your matrix D with a little bit of rearrangement and successively update it.


mrsh (27)
Thank you!i already did it, but now i'm stuck at the sum of the series which is i start with factorial 2!

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
 
#include <fstream.>
#include <iostream>
#include <conio.h>
#include <iomanip>
#include <math.h>
#define n 2
//#define p 2n+1

using namespace std;

void main()
{

	double A[n + 1][n + 1], z[n + 1][n + 1], D[n + 1][n + 1], X[n + 1][n + 1],F[n+1][n+1],factorial = 1;
	int i, j, k, c, d, h = 2 * n + 1, q = 1,sum=0;
	cout.setf(ios::fixed);
	cout.precision(4);

	A[1][1] = 1.0000; A[1][2] = 0.5000;   //A[1][3] = 0.1250;
	A[2][1] = 0.0000; A[2][2] = 0.8750;   //A[2][3] = 0.2500;
	///A[3][1] = 0.0000; //A[3][2] = 0.0000;   //A[3][3] = 0.1250;

	//=================matrix order 1===================

	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= n; j++)
		{
			D[i][j] = A[i][j];
		}
	}

	//==============calculation for matrix multiplication======

	for (k = 2; k <= h; k++)
	{
		factorial = factorial*k;

		for (i = 1; i <= n; i++)
		{
			for (j = 1; j <= n; j++)
			{
				z[i][j] = 0;
				X[i][j] = 0;
				F[i][j] = 0;

				for (int c = 1; c <= n; c++)
				{

					z[i][j] += D[i][c] * A[c][j];
				}

				D[i][j] = z[i][j];
				X[i][j] = D[i][j] / factorial;
				//F[i][j] = X[i][j];
				cout << "   " << X[i][j];

				

				F[i][j]= F[i][j] + X[i][j];
				cout << "   " << F[i][j];

			
	
			}cout << endl;
			
		}cout << endl;

}


	_getch();

}
lastchance (6980)
@mrsh

Presumably your sum is in F[][]. You keep resetting this to zero on line 45. It should be initialised before any of your other matrix loops. If you want an exponential matrix in the end then it should be initialised to the identity matrix.

You are also not doing what was suggested in my previous post. You do not need to store factorial and you could easily avoid having both D and X matrices.
Topic archived. No new replies allowed.