本程序使用java语言在eclipse环境中完成,

一.矩阵存储

此问题中矩阵为501*501的带状矩阵,有5条带,所以用a[5][501]的数组存储.原矩阵中元素a_i,j对应新矩阵中的a_{i-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;    }

五.原点平移原理

幂法:求绝对值最大的特征根.

反幂法:求绝对值最小的特征根.

原点平移原理在本题中大量使用.由题干知,

λ₁<λ₂<λ₃<λ₄<……<λ₅₀₀<λ₅₀₁

通过幂法求得λ_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)₂=λ_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.700113615021694lambdaMin=-0.005557910794214687lambda=9.724634099220056-10.700113615021694 9.7246340992200561 -10.1829340331461752 -9.585707425067623 -9.1726724239280294 -8.6522840078975545 -8.0934838086753786 -7.6594054076924217 -7.1196846486911658 -6.6117643393973239 -6.06610322659509210 -5.58510105262838411 -5.11408352981220412 -4.57887217686512613 -4.09647092625967114 -3.554211215750797715 -3.041090018133268316 -2.53397031113018417 -2.00323076956350918 -1.503557611227384719 -0.993558606007544820 -0.487042673884957121 0.02231736249574780422 0.53241747420686323 1.052898962693453524 1.589445881880881425 2.060330460274292526 2.558075597072831327 3.080240509307067528 3.613620867692304329 4.091378510450617530 4.60303537827914431 5.13292428389843532 5.59490634808324233 6.08093385702686934 6.680354092111584535 7.29387744812679236 7.717111714235640537 8.225220014050238 8.6486660651934839 9.254200344575codition=1925.204273908031det=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也可以进行数学编辑.

仿宋是一种类似瘦金体的字体,非常漂亮,一用仿宋顿时显得高大上.

 

我还是老想记录自己的一点一滴.