Расчет инверсной матрицы с использованием итерационного метода

Question

Мой профессор назначил нам лабораторию, в которой я был в тупик, и не могу найти ничего для помощи в Интернете. Это специально для матрицы с собственным значением < 1.

Учитывая матрицу 3x3 A = {(.5, -1, 0), (0, .6666, 0), (.5, -1, .6666))}, вы должны вычислить обратный, используя A^-1 = I + B + B^2 + B^3 ..... B^20, где I - единичная матрица и B = I - A. Есть два выхода, один из которых A^-1, а другой - A * A^-1 для проверки результата.

До сих пор я пытался использовать две разные функции для умножения. Извините за мое страшное имя. i, j и k - счетчик строк, столбцов и циклов соответственно.

void mPower(double x[][3], double y[][3], double z[][3]) 
{ 
    int i, j, k; 

    for (i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) { 
     for (k = 0; k < 3; k++) 
     z[i][j] = z[i][j] + (x[i][k] * y[k][j]); 
} 
} 

mPower - моя, вероятно, неправильная функция для умножения. В качестве альтернативы я пробовал

void multiplyMatrix(double x[][3], double y[][3], double z[][3]) 
{ 
    int i, j, k; 
    int sum = 0; 

for (i = 0; i <= 2; i++) { 
    for (j = 0; j <= 2; j++) { 
    sum = 0; 
    for (k = 0; k <= 2; k++) { 
     sum = sum + x[i][k] * y[k][j]; 
    } 
    z[i][j] = sum; 
    } 
    } 

} 

Однако я получаю только нули.

void printIt(double a[][3], double b[][3]) 
{ 
int i, j; 

for (i = 0; i < 3; i++) 
{ 
    printf ("\n\t\t\t"); 
    for (j = 0; j < 3; j++) 
     printf ("%7.4f", a[i][j]); 
    printf ("\t"); 
    for (j = 0; j < 3; j++) 
     printf ("%7.4f", b[i][j]); 
     printf ("\n\n\t"); 
} 
} 

int main() 
{ 
    int i, j, k; 
    double b1[][3] = {.5, 1., 0., 0., .6666, 0., -.5, -1., .6666}; 
    double b3[3][3] = {0}; 
    double addOn[3][3]; 
    double aInverse[3][3]; 
    double identity[][3] = {1., 0., 0., 0., 1., 0., 0., 0., 1.}; 
    mPower (b1, b1, b3); 
    for (k = 0; k < 20; k++)  { 
     addOn[i][j] = addOn[i][j] + b3[i][j]; 
     mPower (b1, addOn, b3); 
} 

for (i = 0; i < 3; i++) { 
    for (j = 0; j < 3; j++) 
     aInverse[i][j] = identity[i][j] + addOn[i][j]; 
} 

printIt(b3, aInverse); 
system("pause"); 

return 0; 
} 

PrintIt имеет по крайней мере, сделал форматирование дисплея, чтобы быть хорошо, но я беспокоюсь о для петель в моей главной проблеме. Я просто не могу окунуться в голову, что делать. Я получаю ценности повсюду, и мне действительно трудно понять математику, участвующую в такой проблеме ниши. Любая помощь будет оценена по достоинству.

Answer 1

Я попытался исправить вашу логику (и полученный код очень неэффективен и может быть намного лучше написан). Некоторые вопросы, которые я думаю, что я нашел являются:

1. Переменная аддон не была инициализирована
2. Ваш б было неправильно (плохо вычитание) дано определение B = I-А в своем предисловии.
3. Ваше накопление addOn = addOn + b3 не было в двойном цикле, используя неинициализированные i и j.
4. Вы, кажется, не имеют обратной проверки (а * а^-1) = I

(Обратите внимание, я получаю -0.000, потому что мы используем IEEE двойной точности так, математика не является точным . особенно, когда дело доходит до трети)

#include <stdio.h> 

void mPower(double x[3][3], double y[3][3], double z[3][3]) 
{ 
    int i, j, k; 

    for (i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) { 
     z[i][j] = 0; 
     for (k = 0; k < 3; k++) 
     z[i][j] += (x[i][k] * y[k][j]); 
} 
} 

void mAdd(double x[3][3], double y[3][3]) 
{ 
    int i, j; 
    for (i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) 
     x[i][j] += y[i][j]; 
} 

void mCopy(double x[3][3], double y[3][3]) 
{ 
    int i, j; 
    for (i = 0; i < 3; i++) 
    for (j = 0; j < 3; j++) 
     x[i][j] = y[i][j]; 
} 

void printIt(double a[3][3], double b[3][3]) 
{ 
int i, j; 

for (i = 0; i < 3; i++) 
{ 
    printf ("\n\t\t\t"); 
    for (j = 0; j < 3; j++) 
     printf ("%7.4f", a[i][j]); 
    printf ("\t"); 
    for (j = 0; j < 3; j++) 
     printf ("%7.4f", b[i][j]); 
     printf ("\n\n\t"); 
} 
} 

int main() 
{ 
    int i, j, k; 
    double a[3][3] = {.5, -1., 0., 0., .6666, 0., .5, -1., .6666}; 
    double b1[3][3] = {.5, 1., 0., 0., 1-.6666, 0., -.5, 1., 1-.6666}; 
    double b3[3][3] = {0}; 
    double addOn[3][3] = {0}; 
    double aInverse[3][3] = {0}; 
    double identity[][3] = {1., 0., 0., 0., 1., 0., 0., 0., 1.}; 
printIt(aInverse, b1); 
    mAdd(aInverse, identity); 
printIt(aInverse, b1); 
    mCopy(b3, b1); 
    mCopy(addOn, b1); 
    for (k = 0; k < 20; k++)  { 
    mAdd(aInverse, addOn); 
    mPower(b3, b1, addOn); 
    mCopy(b3, addOn); 
printf("loop iteration %d\n", k); 
printIt(aInverse, addOn); 
} 
printf("aInverse\n"); 
printIt(b1, aInverse); 
printf("verify\n"); 
mPower(a,aInverse,b3); 
printIt(a,b3); 

return 0; 
} 

Результаты я получил:

aInverse 

     0.5000 1.0000 0.0000 2.0000 3.0003 0.0000 


     0.0000 0.3334 0.0000 0.0000 1.5002 0.0000 


     -0.5000 1.0000 0.3334 -1.5001 0.0000 1.5002 

verify 

     0.5000-1.0000 0.0000 1.0000-0.0000 0.0000 


     0.0000 0.6666 0.0000 0.0000 1.0000 0.0000 


     0.5000-1.0000 0.6666 0.0000 0.0000 1.0000