Экспонента матрицы Опции

[Ozzя]

Нашел формулу.

Сообщение отредактировано: Ozzя - 6.12.2005 15:59

[Ozzя]

Для матрицы

A =

    1    1    1
    1    1    1
    1    1    1

Точность невелика, результаты, сопоставимые с матлаб, достигают только при кол-ве итераций, равных 13.

Результаты, полученные программой:

7.3620          6.3620          6.3620
6.3620          7.3620          6.3620
6.3620          6.3620          7.3620

Результаты, полученные Matlab:

  7.3618    6.3618    6.3618
    6.3618    7.3618    6.3618
    6.3618    6.3618    7.3618

Кстати, volvo, в процедуре возведения матрицы в степень - ошибка. Должно быть
For i := 1 To pow-1 Do

{ Порядок матрицы }
Uses Crt;
Const
  size = 3;
  eps = 10E-6;

Type
  TMatrix = Array [1 .. size, 1 .. size] Of real;

Function Factorial(I:LongInt):LongInt;
var
  j,fact:LongInt;
begin
  if i=1 then
    fact:=1
  else
    begin
     fact:=1;
     for j:=2 to i do fact:=fact*j;
   end;
Factorial:=fact;
end;

(* Умножение матриц *)
Procedure matrixMult(Var m: TMatrix; a, b: TMatrix);
Var
  i, j, k: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    Begin
      m[i, j] := 0;
        For k := 1 To size Do
          m[i, j] := m[i, j] + a[i, k] * b[k, j]
    End;
End;

(* Возведение матрицы в степень *)
Procedure matrixPower(Var m: TMatrix; a: TMatrix; pow: Integer);
Var
  i, j: Integer;
  T: TMatrix;
Begin
  If pow = 0 Then
    Begin
      For i := 1 To size Do
      For j := 1 To size Do
        m[i, j] := Byte(i = j);
      Exit
    End;
    move(a, T, SizeOf(T));
    For i := 1 To pow-1 Do
      matrixMult(T, T, a);
      move(T, m, SizeOf(T))
End;

(* Сложение матриц *)
Procedure matrixAdd(Var m: TMatrix; a, b: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    m[i, j] := a[i, j] + b[i, j]
End;

(* Умножение матрицы на число *)
Procedure matrixScale(Var m: TMatrix; a: TMatrix; f: real);
Var
  i, j, k: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
    m[i, j] := a[i, j] / f
End;

(* Печать матрицы *)
Procedure matrixPrint(a: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
    Begin
      For j := 1 To size Do
        write(a[i, j]:16:4);
      WriteLn
    End
End;

Procedure SingleMatrix(Var a: TMatrix);
Var
  i, j: Integer;
Begin
  For i := 1 To size Do
  For j := 1 To size Do
  if i=j then
    a[i, j]:=1
  else
   a[i, j]:=0;
End;

Var
  a,EMatrix,a2: TMatrix;
  T2, Res: TMatrix;
  i, j,n: Integer;
  ExpMatr:TMatrix;
  fact2:real;
begin
  clrscr;
  randomize;
  For i := 1 To size Do
  For j := 1 To size Do
    begin
    a2[i, j] :=2;
    a[i,j]:=a2[i,j];
    end;
  WriteLn('Исходная матрица');
  matrixPrint(a2);

  { Single E-Matrix here }
  SingleMatrix(EMatrix);
{  MatrixAdd(ExpMatr,EMatrix,ExpMatr);
  MatrixAdd(ExpMatr,A,ExpMatr);}

  n:=2;
  repeat
    For i := 1 To size Do
      For j := 1 To size Do
         a2[i,j]:=a[i,j];


    writeln('Matrica v stepeni ',n);
    MatrixPower(a2, a2, n);
    matrixPrint(a2);
    writeln;
    writeln('Matrica v stepeni,podelennay na fact');
    fact2:=Factorial(n);
    MatrixScale(a2,a2,Fact2);
    matrixPrint(a2);
    writeln;
    writeln('Matrica v stepeni podelennay na fact + suuma');
    MatrixAdd(ExpMatr,ExpMatr,A2);
    matrixPrint(ExpMatr);
    writeln;
    n:=n+1;
{    readkey;}
  until n>13;
  MatrixAdd(ExpMatr,EMatrix,ExpMatr);
  MatrixAdd(ExpMatr,A,ExpMatr);
  {abs(delta) <= eps;}

  writeln;
  clrscr;
  matrixPrint(ExpMatr);
  WriteLn('Factorial from ',n,' = ',Factorial(n));
  readln;
end.