本程序使用java语言在eclipse环境中完成,
一.矩阵存储
此问题中矩阵为501*501的带状矩阵,有5条带,所以用a[5][501]的数组存储.原矩阵中元素ai,j对应新矩阵中的ai-j+upCount,j,式中upCount是指对角线上方有多少条带.
为了描述带状矩阵,本程序中定义了BandMatrix类,此类专门处理带状矩阵相关问题.
BandMatrix为带状矩阵,leftBandCount表示左半部分有几条带,upBandCount表示上半部分有几条带,n表示此矩阵为n*n的方阵,a[][]为double类型的数组,get(i,j)函数是获得第i行,第j列函数.set(i,j,x)用于设置第i行,第j列的元素值.mul(Vector)为向量乘法(也就是内积运算).getMaxRoot用幂法获取最大特征根,getMinRoot用反幂法获取最小特征根,getDet()用杜力特尔分解法求行列式值.dolittle为带状矩阵LU分解,分解所得LU存储在同一个矩阵里面,LU矩阵也是带状矩阵.solveEquationByDolittle()函数用于求解线性方程组,它的第一个参数BandMatrix为LU矩阵,Vector为b向量,此函数返回未知数向量.subE(double)函数作用是返回A-xE,返回矩阵也是带状矩阵.
二.向量
本程序定义Vector类来处理向量相关问题.n表示向量的长度,u是一个double数组,getRandomVector(n)静态方法获得一个随机向量,get2Norm()获得向量的二范数,也就是向量的模长.toOne()函数用于返回归一化的向量,mul(Vector)用于向量乘法,返回两个向量的内积.toString()将此Vector转化为String,便于调试.
二.矩阵初始化
void init(BandMatrix matrix) { double[][] a = matrix.a; double b = 0.16, c = -0.064; for (int i = 1; i <= matrix.n; i++) { a[2][i - 1] = (1.64 - 0.024 * i) * sin(0.2 * i) - 0.64 * exp(0.1 / i); } for (int i = 0; i < 501 - 2; i++) { a[0][i + 2] = a[4][i] = c; } for (int i = 0; i < 501 - 1; i++) { a[1][i + 1] = a[3][i] = b; } }
三.幂法求解绝对值最大特征根
// 幂法求解绝对值最大特征根 double getMaxRoot() { double beta = 0; Vector u = Vector.getRandomVector(n); BandMatrix A = this; while (true) { Vector y = u.toOne(); u = A.mul(y); double beta2 = y.mul(u); if (abs(beta - beta2) < Main.epsilon) break; beta = beta2; } return beta; }
四.反幂法求解绝对值最小特征根
// 反幂法求解绝对值最小特征根 double getMinRoot() { double beta = 0; Vector u = Vector.getRandomVector(n); BandMatrix lu = dolittle(); while (true) { Vector y = u.toOne(); u = BandMatrix.solveEquationByDolittle(lu, y); double beta2 = y.mul(u); if (abs(beta - beta2) < Main.epsilon) break; beta = beta2; } return 1 / beta; }
五.原点平移原理
幂法:求绝对值最大的特征根.
反幂法:求绝对值最小的特征根.
原点平移原理在本题中大量使用.由题干知,
λ1<λ2<λ3<λ4<……<λ500<λ501
通过幂法求得λmax,如果λmax<0,说明λmax=λ1,否则说明λmax=λ501.
通过求B=A-λmaxE矩阵的λ’max,则将离λmax最远的那个特征根求出来了.如果λmax=λ1,则λ’max=λ501. 如果λmax=λ501,则λ’max=λ1.
求A与数x最接近的特征值,也需要用到原点平移原理.求矩阵B=A-xE的最小特征根,则将离x最近的特征根找到了.
六.条件数
矩阵的条件数定义为
cond(A)2=λmax/λmin
上式中, λmax和λmin分别表示绝对值最大的特征根和绝对值最小的特征根.
七.求行列式
求矩阵的行列式可以使用高斯消元法,也可以使用杜力特尔分解法.本程序中使用杜力特尔分解法求矩阵行列式值.
八.程序主体
public Main() { BandMatrix A = new BandMatrix(2, 2, col); init(A); double lambdaMax = A.getMaxRoot(); double lambdaMin = A.getMinRoot(); System.out.println("lambdaMax=" + lambdaMax); System.out.println("lambdaMin=" + lambdaMin); BandMatrix B = A.subE(lambdaMax); double lambda = B.getMaxRoot() + lambdaMax; System.out.println("lambda=" + lambda); double lambda501, lambda1; if (lambdaMax > 0) { lambda501 = lambdaMax; lambda1 = lambda; } else { lambda501 = lambda; lambda1 = lambdaMax; } System.out.println(lambda1 + " " + lambda501); for (int k = 1; k <= 39; k++) { double uk = lambda1 + k * (lambda501 - lambda1) / 40.0; BandMatrix matrix = A.subE(uk); double la = matrix.getMinRoot()+uk; System.out.println(k + " " + la); } System.out.println("codition=" + lambdaMax / lambdaMin); System.out.println("det=" + A.getDet()); }
九.程序运行结果
lambdaMax=-10.700113615021694 lambdaMin=-0.005557910794214687 lambda=9.724634099220056 -10.700113615021694 9.724634099220056 1 -10.182934033146175 2 -9.58570742506762 3 -9.172672423928029 4 -8.652284007897554 5 -8.093483808675378 6 -7.659405407692421 7 -7.119684648691165 8 -6.611764339397323 9 -6.066103226595092 10 -5.585101052628384 11 -5.114083529812204 12 -4.578872176865126 13 -4.096470926259671 14 -3.5542112157507977 15 -3.0410900181332683 16 -2.533970311130184 17 -2.003230769563509 18 -1.5035576112273847 19 -0.9935586060075448 20 -0.4870426738849571 21 0.022317362495747804 22 0.532417474206863 23 1.0528989626934535 24 1.5894458818808814 25 2.0603304602742925 26 2.5580755970728313 27 3.0802405093070675 28 3.6136208676923043 29 4.0913785104506175 30 4.603035378279144 31 5.132924283898435 32 5.594906348083242 33 6.080933857026869 34 6.6803540921115845 35 7.293877448126792 36 7.7171117142356405 37 8.2252200140502 38 8.64866606519348 39 9.254200344575 codition=1925.204273908031 det=2.77278614175212E118
十.程序源代码
import static java.lang.Math.abs; import static java.lang.Math.exp; import static java.lang.Math.max; import static java.lang.Math.min; import static java.lang.Math.random; import static java.lang.Math.sin; import static java.lang.Math.sqrt; class Vector { int n; double[] u; public Vector(int n) { this.n = n; u = new double[n]; } // 获取一个随机向量 static Vector getRandomVector(int n) { Vector vector = new Vector(n); for (int i = 0; i < n; i++) { vector.u[i] = random(); } return vector; } // 获取2范式 double get2Norm() { return sqrt(mul(this)); } // 归一化操作 Vector toOne() { Vector v = new Vector(n); double norm = get2Norm(); for (int i = 0; i < u.length; i++) { v.u[i] = u[i] / norm; } return v; } // 向量乘法 double mul(Vector v) { assert (v.n == n); double ans = 0; for (int i = 0; i < u.length; i++) { ans += u[i] * v.u[i]; } return ans; } @Override public String toString() { StringBuilder builder = new StringBuilder("[" + u[0]); for (int i = 1; i < u.length; i++) { builder.append(',').append(u[i]); } builder.append("]"); return builder.toString(); } } // 带状矩阵 class BandMatrix { // leftBandCount表示左半部分有几条带,upBandCount表示右半部分有几条带,n表示方阵的行数和列数 int leftBandCount, upBandCount, n; double[][] a; public BandMatrix(int leftBandCount, int upBandCount, int n) { a = new double[leftBandCount + upBandCount + 1][n]; this.leftBandCount = leftBandCount; this.upBandCount = upBandCount; this.n = n; } double get(int i, int j) { if (i >= j && i - j <= leftBandCount || i < j && j - i <= upBandCount) return a[i - j + upBandCount][j]; else return 0; } void set(int i, int j, double x) { a[i - j + upBandCount][j] = x; } Vector mul(Vector v) { assert (n == v.n); Vector ans = new Vector(v.n); for (int i = 0; i < n; i++) { for (int j = max(0, i - leftBandCount); j <= min(i + upBandCount, n - 1); j++) { ans.u[i] += get(i, j) * v.u[j]; } } return ans; } // 幂法求解绝对值最大特征根 double getMaxRoot() { double beta = 0; Vector u = Vector.getRandomVector(n); BandMatrix A = this; while (true) { Vector y = u.toOne(); u = A.mul(y); double beta2 = y.mul(u); if (abs(beta - beta2) < Main.epsilon) break; beta = beta2; } return beta; } // 反幂法求解绝对值最小特征根 double getMinRoot() { double beta = 0; Vector u = Vector.getRandomVector(n); BandMatrix lu = dolittle(); while (true) { Vector y = u.toOne(); u = BandMatrix.solveEquationByDolittle(lu, y); double beta2 = y.mul(u); if (abs(beta - beta2) < Main.epsilon) break; beta = beta2; } return 1 / beta; } double getDet() { BandMatrix matrix = dolittle(); double ans = 1; for (int i = 0; i < n; i++) { ans *= matrix.get(i, i); } return ans; } BandMatrix dolittle() { BandMatrix lu = new BandMatrix(leftBandCount, upBandCount, n); for (int i = 0; i < n; i++) { for (int j = max(0, i - leftBandCount); j <= i; j++) { double s = get(i, j); for (int k = max(0, i - leftBandCount); k < j; k++) { s -= lu.get(i, k) * lu.get(k, j); } lu.set(i, j, s); } for (int j = i + 1; j <= min(n - 1, i + upBandCount); j++) { double s = get(i, j); for (int k = max(0, i - leftBandCount); k < i; k++) { s -= lu.get(i, k) * lu.get(k, j); } lu.set(i, j, s / lu.get(i, i)); } } return lu; } static Vector solveEquationByDolittle(BandMatrix lu, Vector b) { assert (lu.n == b.n); Vector y = new Vector(b.n); for (int i = 0; i < lu.n; i++) { y.u[i] = b.u[i]; for (int j = max(0, i - lu.leftBandCount); j < i; j++) { y.u[i] -= lu.get(i, j) * y.u[j]; } y.u[i] /= lu.get(i, i); } for (int i = lu.n - 1; i >= 0; i--) { for (int j = i + 1; j <= min(lu.n - 1, i + lu.upBandCount); j++) { y.u[i] -= lu.get(i, j) * y.u[j]; } } return y; } // A-xE BandMatrix subE(double x) { BandMatrix ans = new BandMatrix(leftBandCount, upBandCount, n); for (int i = 0; i < leftBandCount + upBandCount + 1; i++) { for (int j = 0; j < n; j++) { ans.a[i][j] = a[i][j]; } } for (int i = 0; i < n; i++) { ans.a[upBandCount][i] -= x; } return ans; } } public class Main { final static int col = 501; final static double epsilon = 10e-12; void init(BandMatrix matrix) { double[][] a = matrix.a; double b = 0.16, c = -0.064; for (int i = 1; i <= matrix.n; i++) { a[2][i - 1] = (1.64 - 0.024 * i) * sin(0.2 * i) - 0.64 * exp(0.1 / i); } for (int i = 0; i < 501 - 2; i++) { a[0][i + 2] = a[4][i] = c; } for (int i = 0; i < 501 - 1; i++) { a[1][i + 1] = a[3][i] = b; } } void show(BandMatrix matrix) { for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { System.out.print(matrix.get(i, j) + " "); } System.out.println(); } } public Main() { BandMatrix A = new BandMatrix(2, 2, col); init(A); double lambdaMax = A.getMaxRoot(); double lambdaMin = A.getMinRoot(); System.out.println("lambdaMax=" + lambdaMax); System.out.println("lambdaMin=" + lambdaMin); BandMatrix B = A.subE(lambdaMax); double lambda = B.getMaxRoot() + lambdaMax; System.out.println("lambda=" + lambda); double lambda501, lambda1; if (lambdaMax > 0) { lambda501 = lambdaMax; lambda1 = lambda; } else { lambda501 = lambda; lambda1 = lambdaMax; } System.out.println(lambda1 + " " + lambda501); for (int k = 1; k <= 39; k++) { double uk = lambda1 + k * (lambda501 - lambda1) / 40.0; BandMatrix matrix = A.subE(uk); double la = matrix.getMinRoot()+uk; System.out.println(k + " " + la); } System.out.println("codition=" + lambdaMax / lambdaMin); System.out.println("det=" + A.getDet()); } public static void main(String[] args) { new Main(); } }
十一.我说
还是不会使用数学编辑,有时间一定要学会怎样编辑数学公式,照着书上编辑公式,练习完一本书之后肯定就会了,学学tex和word.markdown也可以进行数学编辑.
仿宋是一种类似瘦金体的字体,非常漂亮,一用仿宋顿时显得高大上.
我还是老想记录自己的一点一滴.
441
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
