1.c

来自「实现 了 数值分析 课程 中得 几个 算法」· C语言 代码 · 共 79 行

#include <stdio.h>
#include <math.h>
#define MAX_N 20
#define MAXREPT 100
#define epsilon 0.0001
int main()
{
  int n;
  int i, j, k;
  double err;
  static double a[MAX_N][MAX_N], b[MAX_N][MAX_N], c[MAX_N], g[MAX_N];
  static double x[MAX_N], nx[MAX_N];
  printf("\nInput n value(dim of Ax=c):");
  scanf("%d", &n);
  if (n > MAX_N)
  {
    printf("The input n is larger then MAX_N,please redefine the MAX_N.\n");
    return 1;
  }
  if (n <= 0)
  {
    printf("Please input a number between 1 and %d.\n", MAX_N);
    return 1;
  }
  printf("Now input the matrix a(i,j),i,j=0,...,%d:\n", n - 1);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf", &a[i][j]);

  printf("Now input the matrix b(i),i=0,...,%d:\n", n - 1);
  for (i = 0; i < n; i++)
    scanf("%lf", &c[i]);
  for (i = 0; i < n; i++)
  for (j = 0; j < n; j++)
  {
    b[i][j] =  - a[i][j] / a[i][i];
    g[i] = c[i] / a[i][i];
  }

  for (i = 0; i < MAXREPT; i++)
  {
    for (j = 0; j < n; j++)
      nx[j] = g[j];
    for (j = 0; j < n; j++)
    {
      for (k = 0; k < n; k++)
      {
        if (j == k)
          continue;
        if (j < k)
        {
          nx[j] += b[j][k] *nx[k];
        }
        else
        {
          nx[j] += b[j][k] *x[k];
        }
      }
    }
    err = 0;
    for (j = 0; j < n; j++)
      if (err < fabs(nx[j] - x[j]))
        err = fabs(nx[j] - x[j]);
    for (j = 0; j < n; j++)
      x[j] = nx[j];
    if (err < epsilon)
    {
      printf("Solve...x_i=\n");
      for (i = 0; i < n; i++)
        printf("%f\n", x[i]);
      system("pause");
      return 0;
    }
  }
  printf("After %d repeat,no result...\n", MAXREPT);
  return 1;
}