首页 > 解决方案 > 用于计算矩阵指数的意外输出

问题描述

已编辑

我正在尝试计算矩阵的指数或矩阵的 e。使用一些预期的结果,例如对于 1 + I 的简单复数。我已经实现了一些更改以除以 n 阶乘。仍然有一个类似的问题,即经过一些迭代后,值会像疯了一样爆炸。

输出 这个程序将计算

> 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();
} 

标签: c++classmatrixexp

解决方案


推荐阅读