LU分解法求方程组

LU分解法：当系数矩阵A满足顺序主子式不为0时，可将A分解为为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且分解唯一，然后方程式变为Ly=b,Ux=y,接着先求y，再求出x。
LU公式：

#include "stdafx.h"
#include<stdio.h>
#include<math.h>
int main()
{
  double p,A[12][12],x[12],L[12][12],R[12][12],b[12],y[12];
  int n,i,j,k,m;
  //首先确定阶数（升级代码后）
  printf("请输入方程组阶数(不超过10)：\n");
  scanf("%d",&n);
  //按行输入系数
  for(i=1;i<=n;i++)   //外循环控制行数
  {
    printf("请输入系数矩阵的第%d行：\n",i);
    for(j=1;j<=n;j++)     //内循环控制每行的元素个数
    scanf("%lf",&A[i][j]);
  }
  //因为LU分解法在分解时用不到右端常量，所以单独输入一次就好
  printf("请输入右端向量b:\n");
  for(i=1;i<=n;i++)
      scanf("%lf",&b[i]);
  //根据公式开始进行分解
  for(j=1;j<=n;j++)
      R[1][j]=A[1][j];//（1）上三角矩阵的第一行为A的第一行
  for(i=2;i<=n;i++)
      L[i][1]=A[i][1]/R[1][1];//（2）求出L的第一列
  //内行求解（3）
  for(k=2;k<=n-1;k++)
  {   //k为行
    for(j=k;j<=n;j++)//j=k从主元素开始爬行
    {   p=0.0;
        for(m=1;m<=k-1;m++)
            p=p+L[k][m]*R[m][j]; //确定内行减数 （每次相乘都是一个列加1，一个行加1）
        R[k][j]=A[k][j]-p;    //相减得合阵数
    }
   //内列求解（4）
    for(i=k+1;i<=n;i++)
    {
      p=0.0;
      for(m=1;m<=k-1;m++)
            p=p+L[i][m]*R[m][k]; 
      L[i][k]=(A[i][k]-p)/R[k][k];
    }
  }
  p=0.0;
  for(j=1;j<=n-1;j++)//求R[n][n]
      p=p+L[n][j]*R[j][n];
  R[n][n]=A[n][n]-p;
  y[1]=b[1];             //回代过程求y
  for(k=2;k<=n;k++)
  {
    p=0.0;
    for(j=1;j<=k-1;j++)
    p=p+L[k][j]*y[j];
    y[k]=b[k]-p;
  }
  x[n]=y[n]/R[n][n];//回代过程求方程解x
  for(k=n-1;k>=1;k--)
  {
    p=0.0;
    for(j=k+1;j<=n;j++)
    p=p+R[k][j]*x[j];
    x[k]=(y[k]-p)/R[k][k];
  }
  printf("线性方程组的解为：\n");
  for(i=1;i<=n;i++)
      printf("x[%d]=%.3lf\n",i,x[i]);
  getchar();
  return 0;
}

此方法仍存在大数吃小数