#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;
}
高斯消元和残差计算
发布于 2021-09-21 1 次阅读
Comments NOTHING