數值計算——線性方程組的解法
矩陣分解
矩陣分解是将矩陣拆解為數個矩陣的乘積，可分為三角分解、滿秩分解、QR分解、Jordan分解和SVD(奇異值)分解等，常見的有三種: 1)三角分解法 (LU)：将原正方矩陣分解成一個上三角形矩陣或是排列的上三角形矩陣和一個 下三角形矩陣，用途主要在簡化一個大矩陣的行列式值的計算過程，求 反矩陣，和求解聯立方程組。 2)QR 分解法 (QR) 3)奇異值分解法 (SVD)。
LU分解法的實作
将系數矩陣A轉變成等價兩個矩陣L和U的乘積 ，其中L和U分别是下三角和上三角矩陣。 LU分解在本質上是高斯消元法的一種表達形式。實質上是将A通過初等行變換變成一個上三角矩陣，其變換矩陣就是一個機關下三角矩陣。這正是所謂的杜爾裡特算法:從下至上地對矩陣A做初等行變換，将對角線左下方的元素變成零，然後再證明這些行變換的效果等同于左乘一系列機關下三角矩陣，這一系列機關下三角矩陣的乘積的逆就是L矩陣，它也是一個機關下三角矩陣。
下面是實作算法，其中α和β分别代表L和U中的元素
實作方法： 1）使用matlab：MATLAB以lu函數來執行lu分解法， 其文法為[L,U]=lu(A)。 2）使用java：
package com.kexin.lab2;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;
/**
* LU分解函數進行矩陣分解
*
* @author KeXin
*
*/
public class LU {
/**
* 分解矩陣A = LU
*
* @param a
* @return
*/
@SuppressWarnings({ "rawtypes", "unchecked" })
public static List myLU(double[][] a) {
double[][] U = a; // a是要分解的矩陣
double[][] L = new double[a.length][a.length];
for (int i = 0; i < a.length; i++) {
for (int j = 0; j < a.length; j++) {
L[i][j] = 1;
}
}
for (int j = 0; j < a[0].length - 1; j++) {
if (a[j][j] == 0) {
break;
}
for (int i = j + 1; i < a.length; i++) {
double mult = a[i][j] / a[j][j];
for (int k = j; k < a[i].length; k++) {
U[i][k] = a[i][k] - a[j][k] * mult;
// 得出上三角矩陣U,通過減去矩陣的第一行,第二行,第一行(第二行)得到上三角矩陣
}
L[i][j] = mult; // 得到下三角矩陣是得出上三角矩陣的乘積因子
}
}
List aa = new ArrayList<>();
aa.add(L);
aa.add(U);
return aa;
}
/**
* 回代求解上三角矩陣函數
*
* @param u
* @param b
* @return
*/
public static double[] UpT(double[][] u, double[] b) {
int len = u.length;
double[] x = new double[len];
for (int i = len - 1; i >= 0; i--) {
if (u[i][i] == 0)
break;
x[i] = b[i] / u[i][i];
for (int k = 0; k <= i - 1; k++) {
b[k] = b[k] - u[k][i] * x[i];
}
}
return x;
}
/**
* 回代求解下三角矩陣函數
*
* @param l
* @param b
* @return
*/
public static double[] LowT(double[][] l, double[] b) {
int len = l.length;
double[] x = new double[len];
for (int i = 0; i < len; i++) {
if (l[i][i] == 0)
break;
x[i] = b[i] / l[i][i];
for (int k = 0; k <= i - 1; k++) {
b[k] = b[k] - l[k][i] * x[i];
}
}
return x;
}
/**
* Ax=b 轉化成LUx=b, 先求Ly=b,再求Ux=y
*
* @param l
* @param u
* @param b
* @return
*/
public static double[] myLU(double[][] l, double[][] u, double[] b) {
double[] y = LowT(l, b);
double[] x = UpT(u, y);
return x;
}
@SuppressWarnings("rawtypes")
public static void main(String[] args) {
// 初始化矩陣
double[][] aa = { { 0, 85, 56, 54 }, { 76, 63, 7, 20 }, { 21, 67, 88, 73 }, { 19.3, 43, 30.2, 29.4 } };
List a = myLU(aa);
double[][] L = (double[][]) a.get(0);
double[][] U = (double[][]) a.get(1);
int len = L.length;
// 列印矩陣L、U
System.out.println("LU矩陣分解");
System.out.println("------L is: -------");
for (int i = 0; i < len; i++) {
for (int j = 0; j < i; j++) {
System.out.print(L[i][j] + "\t");
}
System.out.println("1");
}
System.out.println("------U is: -------\n");
for (int i = 0; i < len; i++) {
for (int v = 0; v < i; v++) {
System.out.print("\t");
}
for (int j = i; j < len; j++) {
System.out.print(U[i][j] + "\t");
}
System.out.println();
}
// 使用LU求解Ax=b
double[] b = { 218, 109, 141, 93.7 };
double[] result = myLU(L, U, b);
System.out.println("------X is: -------\n");
for (int i = 0; i < result.length; i++) {
System.out.print(result[i] + "\t");
}
}
}
Cholesky 分解
Cholesky分解法又稱平方根法，是當A為實對稱正定矩陣時，LU三角分解法的變形，即把一個對稱正定的矩陣表示成一個下三角矩陣L和其轉置的乘積。它要求矩陣的所有特征值必須大于零，故分解的下三角的對角元也是大于零的。
使用java實作
package com.kexin.lab2;
public class Cholesky {
/**
* 生成Hilbert矩陣
*
* @param n
* @return
*/
public static double[][] Hilbert(int n) {
double[][] a = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a[i][j] = (double) 1 / (i + j + 1);
}
}
return a;
}
/**
* Cholesky分解
*
* @param a
* @return
*/
public static double[][] Cholesky(double[][] a) {
int n = a.length;
// 楚列斯基分解
for (int k = 0; k < n; k++) {
a[k][k] = Math.sqrt(a[k][k]);
for (int i = k + 1; i < n; i++) {
a[i][k] = a[i][k] / a[k][k];
}
for (int j = k + 1; j < n; j++) {
for (int i = k + 1; i < n; i++) {
a[i][j] = a[i][j] - a[i][k] * a[j][k];
}
}
}
return a;
}
/**
* @param args
*/
public static void main(String[] args) {
// 生成Hilbert矩陣n=5\6\7\8
int n = 5;
double[][] a = Hilbert(n);
double aa[][] = a;
double x[] = new double[n];
double b[] = new double[n];
double bb[] = new double[n];
double y[] = new double[n];
double r[] = new double[n];
double f[] = new double[n];
// 初始化b
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
b[i] = b[i] + a[i][j];
}
}
// 初始化bb
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
bb[i] = bb[i] + aa[i][j];
}
}
// 楚列斯基分解
a = Cholesky(a);
// 前代計算
for (int j = 0; j < n; j++) {
if (a[j][j] != 0) {
x[j] = b[j] / a[j][j];
b[j] = x[j];
}
for (int i = j + 1; i < n; i++) {
b[i] = b[i] - a[i][j] * x[j];
}
}
// 将a進行轉置
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a[i][j] = a[j][i];
}
}
// 回代計算
for (int j = n - 1; j >= 0; j--) {
if (a[j][j] != 0) {
x[j] = b[j] / a[j][j];
}
for (int i = 0; i <= j - 1; i++) {
b[i] = b[i] - a[i][j] * x[j];
}
}
// 輸出結果
System.out.println("------結果為：------");
for (int i = 0; i < n; i++) {
System.out.println(x[i]);
}
// 計算殘差
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
y[i] = y[i] + aa[i][j] * x[j];
}
r[i] = bb[i] - y[i];
}
// 輸出殘差
System.out.println("------殘差為：-----");
for (int i = 0; i < n; i++) {
System.out.println(r[i]);
}
// 誤差的無窮範數
for (int i = 0; i < n; i++) {
y[i] = x[i] - 1;
}
if (y[0] < 0)
y[0] = -y[0];
double fanshu = y[0];
for (int i = 1; i < n; i++) {
if (y[i] < 0)
y[i] = -y[i];
if (y[i] > fanshu)
fanshu = y[i];
}
// 輸出誤差的無窮範數
System.out.println("-------誤差的無窮範數為------");
System.out.println(fanshu);
// 計算條件數
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
f[j] = f[j] + aa[j][i];
}
}
}
}