复数矩阵求逆的 C 语言程序
生活随笔
收集整理的這篇文章主要介紹了
复数矩阵求逆的 C 语言程序
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
關于復數矩陣求逆,如果使用 MATLAB,就非常簡單。我們先用一個 MATLAB 的例子來說明,等會要將 C 語言的程序和 MATLAB 的程序進行對比。
close all; clear all; clc;%定義矩陣a為復數矩陣 a = [[4+2*i,3+1*i,4+3*i,5+5*i];[1+7*i,8+2*i,2+2*i,9+3*i];[4+4*i,5+6*i,1+7*i,7+2*i];[6+1*i,7+8*i,1+4*i,2+5*i]];a_real = real(a); %求矩陣實部 a_imag = imag(a); %求矩陣虛部 a_inv = inv(a); %求矩陣的逆可以看到 a_inv 是這樣的:
下面列出 C 語言的矩陣求逆的代碼。由于VS2012中不能使用 complex.h,我們在這里將復數矩陣的實部和虛部分開定義,數值和上面定義的 a 一樣。
完整代碼及其注釋如下:
#include <stdio.h> #include <stdlib.h> #include <memory.h>#define N 4 #define DEBUG 1 //debug label,0即不打印相關結果,非0打印相關輸出結果//計算實數矩陣的逆 void matrix_inverse_LU(double a[][N],double a_inverse[N][N]) {double l[N][N], u[N][N];double l_inverse[N][N], u_inverse[N][N];//double a_inverse[N][N];int i, j, k;double s, t;memset(l, 0, sizeof(l));memset(u, 0, sizeof(u));memset(l_inverse, 0, sizeof(l_inverse));memset(u_inverse, 0, sizeof(u_inverse));memset(a_inverse, 0, sizeof(u_inverse));for (i = 0; i < N;i++) //計算l矩陣對角線{l[i][i] = 1;}for (i = 0;i < N;i++){for (j = i;j < N;j++){s = 0;for (k = 0;k < i;k++){s += l[i][k] * u[k][j];}u[i][j] = a[i][j] - s; //按行計算u值}for (j = i + 1;j < N;j++){s = 0;for (k = 0; k < i; k++){s += l[j][k] * u[k][i];}l[j][i] = (a[j][i] - s) / u[i][i]; //按列計算l值}}for (i = 0;i < N;i++) //按行序,行內從高到低,計算l的逆矩陣{l_inverse[i][i] = 1;}for (i= 1;i < N;i++){for (j = 0;j < i;j++){s = 0;for (k = 0;k < i;k++){s += l[i][k] * l_inverse[k][j];}l_inverse[i][j] = -s;}}#if DEBUGprintf("test l_inverse:\n");for (i = 0; i < N; i++){for (j = 0; j < N; j++){s = 0;for (k = 0; k < N; k++){s += l[i][k] * l_inverse[k][j];}printf("%f ", s);}putchar('\n');} #endiffor (i = 0;i < N;i++) //按列序,列內按照從下到上,計算u的逆矩陣{u_inverse[i][i] = 1 / u[i][i];}for (i = 1;i < N;i++){for (j = i - 1;j >=0;j--){s = 0;for (k = j + 1;k <= i;k++){s += u[j][k] * u_inverse[k][i];}u_inverse[j][i] = -s / u[j][j];}}#if DEBUGprintf("test u_inverse:\n");for (i = 0;i < N;i++){for (j = 0;j < N;j++){s = 0;for (k = 0;k < N;k++){s += u[i][k] * u_inverse[k][j];}printf("%f ",s);}putchar('\n');} #endiffor (i = 0;i < N;i++) //計算矩陣a的逆矩陣{for (j = 0;j < N;j++){for (k = 0;k < N;k++){a_inverse[i][j] += u_inverse[i][k] * l_inverse[k][j];}}}#if DEBUGprintf("test a:\n");for (i = 0; i < N; i++){for (j = 0; j < N; j++){s = 0;for (k = 0; k < N; k++){s += a[i][k] * a_inverse[k][j];}printf("%f ", s);}putchar('\n');} #endif }//矩陣乘法,由于這里計算的是長度N的方陣,所以實際上ROW,MID,COL的值都是N void MulMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int MID, int COL, double (*arr3)[N]) {double sum=0.0;int i,j,m;for (i = 0; i<ROW; i++){for (j = 0; j<COL; j++){sum = 0.0;for (m = 0; m<MID; m++){sum = sum + arr1[i][m] * arr2[m][j];}arr3[i][j] = sum;}} }//矩陣加法 void PlusMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int COL, double (*arr3)[N]) {int i,j;for(i=0;i<N;i++)//控制行{for(j=0;j<N;j++){arr3[i][j]=arr1[i][j]+arr2[i][j];}} }//矩陣減法 void MinusMatrix(const double (*arr1)[N], const double (*arr2)[N], int ROW, int COL, double (*arr3)[N]) {int i,j;for(i=0;i<N;i++)//控制行{for(j=0;j<N;j++){arr3[i][j]=arr1[i][j]-arr2[i][j];}} }void main() {int i, j, k;double a[N][N]; //N表示矩陣維度,為4double a_real[N][N];double a_imag[N][N];double Ainv[N][N],Binv[N][N],BAinv[N][N],BAinvB[N][N],A_P_BAinvB[N][N],A_P_BAinvB_inv[N][N],AinvB[N][N],AinvB_A_P_BAinvB[N][N],AinvB_A_P_BAinvB_inv[N][N];//將a矩陣的實部和虛部分開定義a_real[0][0] = 4;a_real[0][1] = 3;a_real[0][2] = 4;a_real[0][3] = 5;a_real[1][0] = 1;a_real[1][1] = 8;a_real[1][2] = 2;a_real[1][3] = 9;a_real[2][0] = 4;a_real[2][1] = 5;a_real[2][2] = 1;a_real[2][3] = 7;a_real[3][0] = 6;a_real[3][1] = 7;a_real[3][2] = 1;a_real[3][3] = 2;a_imag[0][0] = 2;a_imag[0][1] = 1;a_imag[0][2] = 3;a_imag[0][3] = 5;a_imag[1][0] = 7;a_imag[1][1] = 2;a_imag[1][2] = 2;a_imag[1][3] = 3;a_imag[2][0] = 4;a_imag[2][1] = 6;a_imag[2][2] = 7;a_imag[2][3] = 2;a_imag[3][0] = 1;a_imag[3][1] = 8;a_imag[3][2] = 4;a_imag[3][3] = 5;//這些計算公式來源于這里:https://wenku.baidu.com/view/2de4c1bc284ac850ad024244.htmlmatrix_inverse_LU(a_real,Ainv);matrix_inverse_LU(a_imag,Binv);MulMatrix(a_imag,Ainv,N,N,N,BAinv);MulMatrix(BAinv,a_imag,N,N,N,BAinvB);PlusMatrix(a_real,BAinvB,N,N,A_P_BAinvB);matrix_inverse_LU(A_P_BAinvB,A_P_BAinvB_inv);MulMatrix(Ainv,a_imag,N,N,N,AinvB);MulMatrix(AinvB,A_P_BAinvB_inv,N,N,N,AinvB_A_P_BAinvB_inv);//最后一步要把AinvB_A_P_BAinvB_inv每個數都求相反數,也是上面鏈接的文獻里說的for (i = 0; i < N; i++){for (j = 0; j < N; j++){AinvB_A_P_BAinvB_inv[i][j] = -AinvB_A_P_BAinvB_inv[i][j];}}//輸出a的逆矩陣的實部矩陣printf("test a_inv_real:\n");for (i = 0; i < N; i++){for (j = 0; j < N; j++){printf("%f ", A_P_BAinvB_inv[i][j]);}printf("\n");}//輸出a的逆矩陣的虛部矩陣printf("test a_inv_imag:\n");for (i = 0; i < N; i++){for (j = 0; j < N; j++){printf("%f ", AinvB_A_P_BAinvB_inv[i][j]);}printf("\n");}//使得窗口不要關閉getchar(); }最后輸出結果如下:
通過對比,可以發現該結果中的實部和虛部與 MATLAB 輸出的結果是一致的。這證明我們的程序是正確的。
總結
以上是生活随笔為你收集整理的复数矩阵求逆的 C 语言程序的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 游黄山(一)
- 下一篇: 相亲软件的心灵测试原理,相亲成功率心理测