使用javascript实现矩阵LU分解

版权声明:本文为博主原创文章,欢迎读者转载学习,转载请注明文章出处。 https://blog.csdn.net/qq_37338983/article/details/79918254

在线性代数中,LU分解是将一个矩阵分解为 L(单位下三角矩阵)和 U(上三角矩阵),可用于求解线性方程组、反矩阵和计算行列式。本文结合LU分解,用javascript实现线性方程组的求解;

假设存在线性系统Ax = b,其中A为n*n矩阵, b为1*n向量,先进行A=LU分解,得到 LUx = b, 再令Y = Ux,得到LY = b,求出Y,再由Y = Ux求解x,具体算法实现如下:

/**
 * 用于实现1-6元线性方程组求解
 * @author: Wujinhua
 * */

function BuildArray () {
    this.Arr = [];
}
BuildArray.prototype.matrix = function(n) {
    for (let i = 0; i < n; i++) {
        let nRows = [];
        for (let j = 0; j < n; j++) {
            nRows.push(0);
        }
        this.Arr.push(nRows);
    }

    return this.Arr;
};
BuildArray.prototype.vector = function (n) {
    for (let i = 0; i < n; i++) {
        this.Arr.push(0);
    }
    return this.Arr;
};

function _argumentsToArray( args ) {
    return [].concat.apply( [], Array.prototype.slice.apply(args) );
}

/**
 * lu分解求解线性方程组
 * @param: n*n矩阵arraysA 和 n维向量arraysB
 * @example:
 * let A = [1,2,3,2,5,2,3,1,5];
 * let b = [4,18,20];
 * let x = lu(A, b);
 */

function lu(arraysA, arraysB) {

    let n = arraysA.length;
    let result1 = _argumentsToArray(arraysA);
    let A = undefined;

    if(n === 2) A = mat2(result1);//2*2矩阵
    if(n === 3) A = mat3(result1);//3*3
    if(n === 4) A = mat4(result1);//4*4
    if(n === 5) A = mat5(result1);//5*5
    if(n === 6) A = mat6(result1);//6*6

    let ll = new BuildArray();//初始化矩阵L
    let l = ll.matrix(n);
    let uu = new BuildArray();//初始化U矩阵
    let u = uu.matrix(n);


    //进行U矩阵的第一行赋值
    for (let i = 0; i < n; i++) {
        u[0][i] = A[0][i];
    }

    //进行L矩阵的第一列赋值
    for (let i = 0; i < n; i++) {
        l[i][0] = A[i][0] / u[0][0];
    }


    //
    for (let i = 1; i < n; i++) {
        for (let j = i; j < n; j++) {
            let sum1 = 0.0;
            for (let k = 0; k < i; k++) {
                sum1 += l[i][k] * u[k][j];
            }
            u[i][j] = A[i][j] - sum1;
        }

        if(i !== n) {
            for (let r = i+1; r < n; r++) {
                let sum2 = 0.0;
                for (let k = 0; k < i; k++) {
                    sum2 += l[r][k] * u[k][i];
                }

                l[r][i] = (A[r][i] - sum2) / u[i][i];
            }
        }
    }

    let y = new BuildArray();
    y.vector(n);
    y[0] = arraysB[0];

    for (let m = 0; m < n; m++) {
        let sum3 = 0.0;
        for (let k = 0; k < m; k++) {
            sum3 += l[m][k] * y[k];
            y[m] = arraysB[m] - sum3;
        }
    }

    let xx = new BuildArray();
    let x = xx.vector(n);
    x[n-1] = y[n-1] / u[n-1][n-1];
    for (let j = n - 2; j >= 0; j--) {
        let sum4 = 0.0;
        for (let k = j + 1; k < n; k++) {
            sum4 += u[j][k] * x[k];
            x[j] = (y[j] - sum4) / u[j][j];
        }
    }

    return x;
}

其中用到构建n*n矩阵的构造函数matn(),其代码如下:

function mat2()
{
    var v = _argumentsToArray( arguments );
    var m = [];
    switch ( v.length ) {
    case 0:
        v[0] = 1;
    case 1:
        m = [
            vec2( v[0],  0.0 ),
            vec2(  0.0, v[0] )
        ];
        break;
    default:
        m.push( vec2(v) );  v.splice( 0, 2 );
        m.push( vec2(v) );
        break;
    }
    m.matrix = true;
    return m;
}
//----------------------------------------------------------------------------
function mat3()
{
    var v = _argumentsToArray( arguments );
    var m = [];
    switch ( v.length ) {
    case 0:
        v[0] = 1;
    case 1:
        m = [
            vec3( v[0],  0.0,  0.0 ),
            vec3(  0.0, v[0],  0.0 ),
            vec3(  0.0,  0.0, v[0] )
        ];
        break;
    default:
        m.push( vec3(v) );  v.splice( 0, 3 );
        m.push( vec3(v) );  v.splice( 0, 3 );
        m.push( vec3(v) );
        break;
    }
    m.matrix = true;
    return m;
}
//----------------------------------------------------------------------------
function mat4()
{
    var v = _argumentsToArray( arguments );
    var m = [];
    switch ( v.length ) {
    case 0:
        v[0] = 1;
    case 1:
        m = [
            vec4( v[0], 0.0,  0.0,   0.0 ),
            vec4( 0.0,  v[0], 0.0,   0.0 ),
            vec4( 0.0,  0.0,  v[0],  0.0 ),
            vec4( 0.0,  0.0,  0.0,  v[0] )
        ];
        break;
    default:
        m.push( vec4(v) );  v.splice( 0, 4 );
        m.push( vec4(v) );  v.splice( 0, 4 );
        m.push( vec4(v) );  v.splice( 0, 4 );
        m.push( vec4(v) );
        break;
    }
    m.matrix = true;
    return m;
}
function mat5()
{
    var v = _argumentsToArray( arguments );
    var m = [];
    switch ( v.length ) {
        case 0:
            v[0] = 1;
        case 1:
            m = [
                vec5( v[0], 0.0,  0.0,   0.0,  0.0 ),
                vec5( 0.0,  v[0], 0.0,   0.0,  0.0 ),
                vec5( 0.0,  0.0,  v[0],  0.0,  0.0 ),
                vec5( 0.0,  0.0,  0.0,   v[0], 0.0 ),
                vec5( 0.0,  0.0,  0.0,   0.0,  v[0] )
            ];
            break;

        default:
            m.push( vec5(v) );  v.splice( 0, 5 );
            m.push( vec5(v) );  v.splice( 0, 5 );
            m.push( vec5(v) );  v.splice( 0, 5 );
            m.push( vec5(v) );  v.splice( 0, 5 );
            m.push( vec5(v) );
            break;
    }
    m.matrix = true;
    return m;
}
function mat6()
{
    var v = _argumentsToArray( arguments );
    var m = [];
    switch ( v.length ) {
        case 0:
            v[0] = 1;
        case 1:
            m = [
                vec6( v[0], 0.0,  0.0,   0.0,  0.0,  0.0 ),
                vec6( 0.0,  v[0], 0.0,   0.0,  0.0,  0.0 ),
                vec6( 0.0,  0.0,  v[0],  0.0,  0.0,  0.0 ),
                vec6( 0.0,  0.0,  0.0,   v[0], 0.0,  0.0 ),
                vec6( 0.0,  0.0,  0.0,   0.0,  v[0], 0.0 ),
                vec6( 0.0,  0.0,  0.0,   0.0,  0.0, v[0] )
            ];
            break;
        default:
            m.push( vec6(v) );  v.splice( 0, 6 );
            m.push( vec6(v) );  v.splice( 0, 6 );
            m.push( vec6(v) );  v.splice( 0, 6 );
            m.push( vec6(v) );  v.splice( 0, 6 );
            m.push( vec6(v) );
            break;
    }
    m.matrix = true;
    return m;
}

向量的构造函数vecn()代码如下:

function vec2()
{
    var result = _argumentsToArray( arguments );
    switch ( result.length ) {
    case 0: result.push( 0.0 );
    case 1: result.push( 0.0 );
    }
    return result.splice( 0, 2 );
}

function vec3()
{
    var result = _argumentsToArray( arguments );
    switch ( result.length ) {
    case 0: result.push( 0.0 );
    case 1: result.push( 0.0 );
    case 2: result.push( 0.0 );
    }
    return result.splice( 0, 3 );
}

function vec4()
{
    var result = _argumentsToArray( arguments );
    switch ( result.length ) {
    case 0: result.push( 0.0 );
    case 1: result.push( 0.0 );
    case 2: result.push( 0.0 );
    case 3: result.push( 1.0 );
    }
    return result.splice( 0, 4 );
}


function vec5()
{
    var result = _argumentsToArray( arguments );
    switch ( result.length ) {
        case 0: result.push( 0.0 );
        case 1: result.push( 0.0 );
        case 2: result.push( 0.0 );
        case 3: result.push( 0.0 );
        case 4: result.push( 0.0 );
    }
    return result.splice( 0, 5 );
}

function vec6()
{
    var result = _argumentsToArray( arguments );
    switch ( result.length ) {
        case 0: result.push( 0.0 );
        case 1: result.push( 0.0 );
        case 2: result.push( 0.0 );
        case 3: result.push( 0.0 );
        case 4: result.push( 0.0 );
        case 5: result.push( 0.0 );
    }

    return result.splice( 0, 6 );
}

用一个五元一次方程组测试如下:

 let A1 = mat5(
        2, 3, -5, 6, -7,
        1, 2, -3, 3, -2,
        3, -7, 2, 3, -4,
        4, 3, -5, 6, -7,
        5, 3, -3, 1, -2
    );

    let b = [10, 15, 13, 22, 33];

    let x = lu(A1, b);

结果:

结果与用MATLAB计算结果一致,证明程序可实现线性方程组求解。

如果需要下载本项目代码,请前往我的Github上下载:https://github.com/Brucewu998/LU-Decomposition

没有更多推荐了,返回首页