⁄(⁄ ⁄•⁄ω⁄•⁄ ⁄)⁄. 菜鸟不吃草!

ヾ(o◕∀◕)ノ ヾ(o◕∀◕)ノ ヾ(o◕∀◕)ノ 亲爱的观众老爷们,有错请指出!...

[计算方法]列主元消去法 LU分解法求线性方程解

/*计算方法实验二
 *WYC
 *计科1604
 */
#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;

const int MAXSIZE = 20;
void columnPivoting(double a[MAXSIZE][MAXSIZE],int n, double b[]);   //声明函数columnPivoting

int main(){
    double matrix[MAXSIZE][MAXSIZE];//系数
    double b[MAXSIZE]; //b
    int sizeM;
    cin>> sizeM;
    for(int i=0;i<sizeM;i++){
        for(int j=0;j<sizeM;j++)
            cin>> matrix[i][j];
        cin>>b[i];
    }
    columnPivoting(matrix,sizeM,b);
    return 0;
}

//列主元消去法
void columnPivoting(double a[MAXSIZE][MAXSIZE], int n, double b[]){
    int i, j, k;
    int col = 0, row = 0;
    for (k = 0; k < n - 1; k++){
        double tmp = 0;
        //找出消元列中最大的那个元素所在的位置
        for (i = k; i < n ; i++)
            if (fabs(a[i][k]) > tmp){
                tmp = fabs(a[i][k]);
                //cout << "MAX " << tmp << endl;//
                row = i;
                col = k;
            }

        //如果该对角线元素是0,同样不能用高斯消元法来求解
        if (a[row][row] == 0){
            cout << "can't solve" << endl;
            return;
        }
        //将找出的行进行交换
        if (k != row){
            for (i = 0; i < n; i++){
                swap(a[row][i], a[k][i]);
                swap(b[k], b[row]);
            }
        }
        //消元过程
        double c[n];
        for (j = k + 1; j < n; j++){
            c[j] = a[j][k] / a[k][k];
            //cout << c[j] << endl;
        }
        for (i = k + 1; i < n; i++){
            for (j = 1; j < n; j++){
                a[i][j] = a[i][j] - c[i] * a[k][j];
            }
            b[i] = b[i] - c[i] * b[k];
        }

    }

    double x[n];
    x[n - 1] = b[n - 1] / a[n - 1][n - 1];
    for (i = n - 2; i >= 0; i--){
        double sum = 0;
        for(j = i + 1; j < n; j++)
            sum += a[i][j] * x[j];
        x[i] = (b[i] - sum)/a[i][i];
    }
    for(int i=0;i<n;i++)
        printf("x[%d] = %.4f\n",i+1,x[i]);  //输出结果
    printf("\n");
}

/*
 4
 1 1 0 3 4
 2 1 -1 1 1
 3 -1 -1 3 -3
 -1 2 3 -1 4
 3
 1 1 -1 1
 1 2 -2 0
 -2 1 1 1
 4
 3 1 0 0 1
 2 3 1 0 0 
 0 2 3 1 1
 0 0 1 3 0

 */
/*计算方法实验二 LU Decomposition
 *WYC
 *计科1604
 */
#include<iostream>
using namespace std;

const int MAXSIZE = 20;
void LUDecomposition(double a[MAXSIZE][MAXSIZE], double b[MAXSIZE],int n);

int main()
{
    double a[MAXSIZE][MAXSIZE];
    double b[MAXSIZE];
    int n;
    cin>> n;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            cin>> a[i][j];
        cin>>b[i];
    }
    LUDecomposition(a, b, n);
    return 0;
}

//矩阵的ALU分解
void LUDecomposition(double a[MAXSIZE][MAXSIZE], double b[MAXSIZE],int n)
{
    double l[MAXSIZE][MAXSIZE] = { 0 };
    double u[MAXSIZE][MAXSIZE] = { 0 };
    int i, r, k;
    //进行U的第一行的赋值
    for (i = 0; i<n; i++)
    {
        u[0][i] = a[0][i];
    }

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

    //计算U的剩下的行数和L的剩下的列数
    for (r = 1; r<n; r++)
    {
        for (i = r; i <n; i++)
        {
            double sum1 = 0;
            for (k = 0; k < r; k++)
            {
                sum1 += l[r][k] * u[k][i];
                //cout << "" << r << "" << sum1 << endl;
            }
            u[r][i] = a[r][i] - sum1;
        }


        if(r!=n)
            for(i=r+1;i<n;i++)
            {
                double sum2 = 0;
                for (k = 0; k<r; k++)
                {
                    sum2 += l[i][k] * u[k][r];
                }
                l[i][r] = (a[i][r] - sum2) / u[r][r];
            }

    }

    double y[MAXSIZE] = { 0 };
    y[0] = b[0];
    for (i = 1; i<n; i++)
    {
        double sum3 = 0;
        for (k = 0; k<i; k++)
            sum3 += l[i][k] * y[k];
        y[i] = b[i] - sum3;
    }

    double x[MAXSIZE] = { 0 };
    x[n - 1] = y[n - 1] / u[n - 1][n - 1];
    for (i = n - 2; i >= 0; i--)
    {
        double sum4 = 0;
        for (k = i + 1; k<n; k++)
            sum4 += u[i][k] * x[k];
        x[i] = (y[i] - sum4) / u[i][i];
    }
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]=" << x[i] << endl;

}


阅读更多
版权声明:本文为博主原创文章,未经博主允许不可转载。 https://blog.csdn.net/orion_pistachio/article/details/79971470
个人分类: 计算方法
上一篇【计算方法】实验一 非线性方程求根数值解法
下一篇【计算方法】雅可比迭代法和高斯-赛德尔迭代法求线性方程根
想对作者说点什么? 我来说一句

高斯列主元消去法原理及代码

2010年04月26日 29KB 下载

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

关闭
关闭