高斯消元和残差计算

发布于 25 天前  42 次阅读


#include <stdio.h>
#include <cmath>
#include <cstdlib>
#include <iostream>
using namespace std;
#define N 100
const float eps = 0.0000001;
double b[N],x[N];
float r[N];
float a[N][N];
double A[N][N];
float ans[N];
int n;
int select(int x) //部分选主元
{
    int f = x;
    for (int i = x + 1; i <= n; i++)
        if (fabs(a[i][x]) > fabs(a[f][x]))
            f = i;
    if (f != x)
    {
        for (int i = x; i <= n + 1; i++)
            swap(a[f][i], a[x][i]);
    }
    return 0;
}
int selectFull(int x) //全选主元
{
    int px = x, py = x;
    for (int i = x; i <= n; i++)
        for (int j = x; j <= n; j++)
            if (fabs(a[i][j]) > fabs(a[px][py]))
            {
                px = i;
                py = j;
            }
    if (px != x || py != x)
    {
        for (int i = x; i <= n + 1; i++)
            swap(a[px][i], a[x][i]);
        for (int i = 1; i <= n; i++)
            swap(a[i][py], a[i][x]);
        //swap(a[py][n + 1], a[x][n + 1]);   
    }
    return 0;
}
int Gauss()
{
    float t;
    for (int i = 1; i <= n; i++)
    {
        selectFull(i);
        if (fabs(a[i][i]) <= eps)
        {
            printf("No Answers\n");
            return 0;
        }
        for (int j = i + 1; j <= n; j++)
        {
            if (fabs(a[j][i]) > eps)
            {
                t = a[i][i] / a[j][i];
                for (int k = i; k <= n + 1; k++)
                    a[j][k] = a[j][k] * t - a[i][k];
            }
        }
    }
    //解出答案
    for (int i = n; i; i--)
    {
        for (int j = i + 1; j <= n; j++)
        {
            a[i][n + 1] -= a[i][j] * a[j][n + 1];
        }
        a[i][n + 1] /= a[i][i];
    }
    return 0;
}
int work() //计算残差
{
    for (int i = 1; i <= n; i++)
    {
        b[i] = A[i][n + 1];
        x[i] = a[i][n + 1];
    }
    for (int i = 1; i <= n; i++)
    {
        double t = 0;
        for (int j = 1; j <= n; j++)
             t += A[i][j] * x[j];
        r[i] = b[i] - t; 
    }
    return 0;
}
int main()
{
    scanf_s("%d", &n);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n + 1; j++)
        {
            cin >> a[i][j];
            A[i][j] = a[i][j];
        }
    Gauss();
    swap(a[1][n + 1], a[4][n + 1]);
    swap(a[3][n + 1], a[2][n + 1]);
    for (int i = 1; i <= n; i++)
        printf("%f\n", a[i][n + 1]);
    work();
    printf("残差为:\n");
    for (int i = 1; i <= n; i++)
        printf("%f\n", r[i]);
    return 0;
}

「雪霁融雾月,冰消凝夜雨」