求逆矩阵 —— LU分解法「建议收藏」

求逆矩阵 —— LU分解法「建议收藏」LU分解:算法步骤:1.将A矩阵分解为L下三角矩阵和U上三角矩阵。step1.L对角线填充为1step2.step3.step4.U是按行迭代计算,L是按列迭代计算,UL交错计算,且U先L一步fork=1tom-1:{}2.分别对L和U求逆,得到Linv和Uinv.step1….

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

LU分解:

A=L*U\Rightarrow A^{-1} = (L*U)^{-1}=L^{-1}*U^{-1}

算法步骤:

1. 将 A 矩阵分解为 L 下三角矩阵和 U 上三角矩阵。

step1. L对角线填充为1

step2. U_{0,j}=A_{0,j}(j=0,...,m-1)

step3. L_{i,0}=\frac{A_{i,0}}{U_{0,0}},(i=1,...,m-1)

step4. U是按行迭代计算,L是按列迭代计算,UL交错计算,且U先L一步

    for  k = 1  to  m-1: {

   U_{k,j}=\frac{A_{k,j}-\sum_{t=0}^{k-1}(L_{i,t}U_{t,j})}{L_{k,k}}=A_{k,j}-\sum_{t=0}^{k-1}(L_{k,t} U_{t,j}),(j=k,...,m-1)

L_{i,k}=\frac{A_{i,k}-\sum_{t=0}^{k-1}(L_{i,t}U_{t,k})}{U_{k,k}},(i=k,...,m-1)

}

2. 分别对 L 和 U 求逆,得到 Linv 和 Uinv.

step1. 列顺序行顺序, for  j = 0  to  m-1, i = j  to  m-1

3. Ainv = Linv * Uinv. 

实现代码(java):

public class TestMain {

	public static void main(String[] args) {
		double[][] A = { { 4, 2, 1, 5 }, { 8, 7, 2, 10 }, { 4, 8, 3, 6 }, { 6, 8, 4, 9 } };
		double[][] U = new double[4][4];
		double[][] L = new double[4][4];
		double[][] Uinv = new double[4][4];
		double[][] Linv = new double[4][4];
		print(A);
		int size = 4;
		// fill 1 at L's diagonal
		for (int i = 0; i < size; i++) {
			L[i][i] = 1.0;
		}
		// calculate the U's first row
		for (int j = 0; j < size; j++) {
			U[0][j] = A[0][j];
		}
		// calculate the L's first column
		for (int i = 1; i < size; i++) {
			L[i][0] = A[i][0] / U[0][0];
		}
		// iterative calculation of other rows and columns
		for (int k = 1; k < size; k++) {
			// the k_th row of U
			for (int j = k; j < size; j++) {
				// sum(L_kt*U_tj),t
				double s = 0.0;
				for (int t = 0; t < k; t++) {
					s += L[k][t] * U[t][j];
				}
				// U_kj = A_kj - s
				U[k][j] = A[k][j] - s;
			}
			// the k_th column of L
			for (int i = k; i < size; i++) {
				// sum(L_it*U_tk),t
				double s = 0.0;
				for (int t = 0; t < k; t++) {
					s += L[i][t] * U[t][k];
				}
				// L_ik = (A_ik - s) / U_kk
				L[i][k] = (A[i][k] - s) / U[k][k];
			}
		}
		
		// calculate L_inv
		for (int j = 0; j < size; j++) {
			for (int i = j; i < size; i++) {
				if (i == j) Linv[i][j] = 1 / L[i][j];
				else if (i < j) Linv[i][j] = 0;
				else {
					double s = 0.0;
					for (int k = j; k < i; k++) {
						s += L[i][k] * Linv[k][j];
					}
					Linv[i][j] = -Linv[j][j] * s;
				}
			}
		}
		
		// calculate U_inv
		for (int j = 0; j < size; j++) {
			for (int i = j; i >= 0; i--) {
				if (i == j) Uinv[i][j] = 1 / U[i][j]; 
				else if (i > j) Uinv[i][j] = 0;
				else {
					double s = 0.0;
					for (int k = i + 1; k <= j; k++) {
						s += U[i][k] * Uinv[k][j];
					}
					Uinv[i][j] = -1 / U[i][i] * s;
				}
			}
		}
		
		double[][] inv = new double[4][4];
		for (int i = 0; i < size; i++) {
			for (int j = 0; j < size; j++) {
				for (int k = 0; k < size; k++) {
					inv[i][j] += Uinv[i][k] * Linv[k][j];
				}
			}
		}
		
		print(L);
		print(U);
		print(Linv);
		print(Uinv);
		print(inv);
	}

	public static void print(double[][] M) {
		for (int i = 0; i < M.length; i++) {
			for (int j = 0; j < M[0].length; j++) {
				System.out.printf("%8.2f ", M[i][j]);
			}
			System.out.println();
		}
		System.out.println();
	}

}

程序输出:

    4.00     2.00     1.00     5.00

    8.00     7.00     2.00    10.00

    4.00     8.00     3.00     6.00

    6.00     8.00     4.00     9.00

 

    1.00     0.00     0.00     0.00

    2.00     1.00     0.00     0.00

    1.00     2.00     1.00     0.00

    1.50     1.67     1.25     1.00

 

    4.00     2.00     1.00     5.00

    0.00     3.00     0.00     0.00

    0.00     0.00     2.00     1.00

    0.00     0.00     0.00     0.25

 

    1.00     0.00     0.00     0.00

   -2.00     1.00     0.00     0.00

    3.00    -2.00     1.00     0.00

   -1.92     0.83    -1.25     1.00

 

    0.25    -0.17    -0.13    -4.50

    0.00     0.33    -0.00    -0.00

    0.00     0.00     0.50    -2.00

    0.00     0.00     0.00     4.00

 

    8.83    -3.67     5.50    -4.50

   -0.67     0.33     0.00     0.00

    5.33    -2.67     3.00    -2.00

   -7.67     3.33    -5.00     4.00

 

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/171656.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号