Matrix Exponential

I define my 2D matrix as vector <vector < double> > MXd2 and i want to solve exponential matrix, can I used this https://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html#aeb6c0eb89cc982629305f6c7e0791caf to solve my problem?

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
#include<iostream>
#include <vector>
#include <conio.h>
#include <iomanip>
#define KMAX 2

using namespace std;

void main()
{


	int i, j, M, L = 0, k = 1, A = 0, R = 0, z = 1, C = 0, c = 1, B = 0, q = 1;
	int Vmax = 3, KM = 5;
	double p = 0.95, a1 = 20.79, a2 = 0.005;


	vector<double> P(KMAX+1);
	vector<double> Q(KMAX+1);
	cout.setf(ios::fixed);
	cout.precision(4);

	for (i = 1; i <=KMAX; i++)
	{
		P[i] = (p*Vmax*i) / (KM + i);
		Q[i] = (1 - p)*Vmax*i / (KM + i);
		cout << "P[" << i << "] = " << P[i] << "\t" << "Q[" << i << "] = " << Q[i] << endl;
	}

	/////////////////////////////////////////////////////////////////////////////////////////
	//matrix
	////////////////////////////////////////////////////////////////////////////////////////

	M = ((KMAX + 1)*(KMAX + 2)) / 2;

	vector <double> MX(M, 0);
	vector < vector <double> >MXd2(M, MX);
	for (i = 0; i < MXd2.size(); i++)
	{
		for (j = 0; j < MXd2[i].size(); j++)
			cout << "   " << MXd2[i][j];
		cout << endl;

	}

	for (i = KMAX; i >= 1; i--)// call P in certain position as formulate in row and column below
	{
		for (j = 1; j <= i; j++)
		{
			MXd2[j - 1 + L][KMAX + k] = P[A+1];
			cout << "MXd2[" << j - 1 + L << "][" << KMAX + k << "]=" << MXd2[j - 1 + L][KMAX + k] << "    " << P[A+1] << endl;
			k++;
		}
		cout << endl;
		L = k + A;
		A++;

	}

	for (i = KMAX; i >= 1; i--)// call Q in certain position as formulate in row and column below
	{
		for (j = 1; j <= i; j++)
		{
			MXd2[j + R][KMAX + z] = Q[C+1];
			cout << "MXd2[" << j + R << "][" << KMAX + z << "] = " << MXd2[j + R][KMAX + z] << "    " << Q[C+1] << endl;
			z++;
		}
		cout << endl;
		R = z + C;
		C++;

	}

	for (i = KMAX + 1; i >= 1; i--)// calculation for diagonal

	{
		for (j = 1; j <= i; j++)
		{
			MXd2[B + c - q][B + c - q] = -(a1*double(j - 1) + a2*double(B)*double(B)) - ((Vmax*double(B)) / (KM + double(B)));
			cout << "MXd2[" << B + c - q << "][" << B + c - q << "] = " << MXd2[B + c - q][B + c - q] << endl;
			c++;

		}
		q++;
		B++;
	}


	cout << "Here is the matrix:" << endl; //double array matrix
	for (i = 0; i < M; i++)
	{
		for (j = 0; j < M; j++)
		{
			cout << setw(7) << MXd2[i][j] << "  ";
		}
		cout << endl;
	}

	_getch();
}
	

Last edited on
mrsh

All your program does is set up a matrix MXd2 (which it does in an extremely complex way!)

I'm sure that you can use the decomposition in your link, provided that you install the package Eigen, include the relevant header(s) in your source code and call the specified routine. You aren't being asked to compute this decomposition from scratch yourself (as far as I can see). At least, there is no documentation for how to do it. Your matrix will have to be in a form compatible with the Matrix class of Eigen - you will have to check the documentation for that.

Your matrix is already upper-triangular. That is quite a useful form in itself (e.g. matrix equations involving it are readily solved by back-substitution), so I'm slightly surprised that you want to decompose it.

Just a bit of advice:
- keep it simple! add documentation and appropriate comments, and try to code in a manner where it is easy to see what is going on.
- it should be
int main()
not
void main()
in C++.
Last edited on
closed account (48T7M4Gy)
Putting aside the the matrix being upper triangular for a second, why not write the solver as an exercise instead of using Eigen or any other package.

Even the least optimal reduction methods for solving a matrix take only about 30 or 40 lines to write and if the solver is operating on a 6x6 matrix which is what this one appears to be, processing time is negligible.
Topic archived. No new replies allowed.