迭代法.cpp

largreance插值 迭代法 迭代改善法

#include<stdio.h>
#include<math.h>
void main()
{
  int i,j,k,n,m,r,h;
  double a[20][20],z[20][20],x[20][20],y[20][20],q[20][20];
  double s[20],c,w,p,t[20];
  printf("输入矩阵阶数input n");
  scanf("%d",&n);
  printf("输入矩阵n行n列input a[][]");
  for(i=1;i<=n;i++)//输入
  {
    for(j=1;j<=n;j++)
    {
     scanf("%lf",&a[i][j]);
    }
  } 
  for(i=1;i<=n;i++)//输入
  {
    for(j=1;j<=n;j++)
    {
     q[i][j]=a[i][j];
    }
  } 
  printf("请输入等式右边的值");
    for(j=1;j<=n;j++)//输入
    {
     scanf("%lf",&y[0][j]);
    }
  for(j=1;j<=n;j++)
  {
	  a[j][n+1]=y[0][j];
  }
for(i=1;i<=n;i++)//输出
  {
    for(j=1;j<=n+1;j++)
    {
     printf("%lf",a[i][j]);
    }
	printf("\n");
  }
//求X的近似值
  for(r=1;r<n;r++)//doolitte法求x的值，选主元
  {
	    for(i=r;i<=n;i++)
		{
			s[i]=a[i][r];
			for(k=1;k<r;k++)
				s[i]=s[i]-a[i][k]*a[k][r];
			a[i][r]=s[i];
		}
	    p=0;
		for(i=r;i<=n;i++)//主元是0退出
		{
			if(fabs(s[i])>p)
				p=s[i];
			m=i;
		}
			if(fabs(p)<=0)
			{
				printf("错误");
				break;
			}
			if(m!=r)//换行
				for(j=1;j<=n+1;j++)
				{
					w=a[r][j];
					a[r][j]=a[m][j];
					a[m][j]=w;
				}
				for(i=r+1;i<=n;i++)
					a[i][r]=a[i][r]/a[r][r];
                for(j=r+1;j<=n+1;j++)
				{
					c=a[r][j];
					for(k=1;k<=r-1;k++)
					    c=c-a[r][k]*a[k][j];
					a[r][j]=c;
				}
	}
	c=a[n][n];//计算分解
	for(k=1;k<=n-1;k++)
		c=c-a[n][k]*a[k][n];
	a[n][n]=c;
    
	c=a[n][n+1];
	for(k=1;k<=n-1;k++)
		c=c-a[n][k]*a[k][n+1];
	a[n][n+1]=c;
	
	if(a[n][n]==0)
	printf("error");
	a[n][n+1]=a[n][n+1]/a[n][n];
    for(i=n-1;i>=1;i--)
	{ 
		for(k=i+1;k<=n;k++)
			a[i][n+1]=a[i][n+1]-a[i][k]*a[k][n+1];
		a[i][n+1]=a[i][n+1]/a[i][i];
	}
   for(i=1;i<=n;i++)
   {
    x[1][i]=a[i][n+1];
   } 
   printf("输出近似值");	
   for(i=1;i<=n;i++)
		printf("%f",x[1][i]);
	printf("\n* * * * * * * * * * * * * * * * * * * * * * \n");
//求值成功，并以输出
	h=1;
   w=1;
   while(w==1)
   {
   for(i=1;i<=n;i++)
	   t[i]=0;
   printf("输出剩余向量:");  
   for(i=1;i<=n;i++)
	   {
		   for(j=1;j<=n;j++)
		  {
			  t[i]=q[i][j]*x[h][j]+t[i];
		  }		  			
		  y[h][i]=y[0][i]-t[i];
	   }
 for(i=1;i<=n;i++)//输出 
		printf("%lf  ",y[h][i]);
 	   printf("\n");
	   for(j=1;j<=n;j++)
  {
	  a[j][n+1]=y[h][j];
  }
	   for(r=1;r<n;r++)//doolitte法选主元
  {//求修正向量Z
	    for(i=r;i<=n;i++)
		{
			s[i]=a[i][r];
			for(k=1;k<r;k++)
				s[i]=s[i]-a[i][k]*a[k][r];
			a[i][r]=s[i];
		}
	    p=0;
		for(i=r;i<=n;i++)//主元是0退出
		{
			if(fabs(s[i])>p)
				p=s[i];
			m=i;
		}
			if(fabs(p)<=0)
			{
				printf("错误");
				break;
			}
			if(m!=r)//换行
				for(j=1;j<=n+1;j++)
				{
					w=a[r][j];
					a[r][j]=a[m][j];
					a[m][j]=w;
				}
				for(i=r+1;i<=n;i++)
					a[i][r]=a[i][r]/a[r][r];
                for(j=r+1;j<=n+1;j++)
				{
					c=a[r][j];
					for(k=1;k<=r-1;k++)
					    c=c-a[r][k]*a[k][j];
					a[r][j]=c;
				}
	}
	c=a[n][n];//计算分解
	for(k=1;k<=n-1;k++)
		c=c-a[n][k]*a[k][n];
	a[n][n]=c;
    
	c=a[n][n+1];
	for(k=1;k<=n-1;k++)
		c=c-a[n][k]*a[k][n+1];
	a[n][n+1]=c;
	
	if(a[n][n]==0)
	printf("error");
	a[n][n+1]=a[n][n+1]/a[n][n];
    for(i=n-1;i>=1;i--)//最后计算 
	{ 
		for(k=i+1;k<=n;k++)
			a[i][n+1]=a[i][n+1]-a[i][k]*a[k][n+1];
		a[i][n+1]=a[i][n+1]/a[i][i];
	}
   for(i=1;i<=n;i++)//z=a[n]
   {
    z[h][i]=a[i][n+1];
   } 
   printf("求修正向量Z:");
   for(i=1;i<=n;i++)//输出 
		printf("%lf ",z[h][i]);
   for(i=1;i<=n;i++)//x=x+z
   {
	   x[h+1][i]=x[h][i]+z[h][i];
   }
	   h++;
	   if(z[h][1]<0.01)
		   w=0;
	   else w=1;
   }
   printf("\n输出改进解:");
   for(i=1;i<=n;i++)//输出 
		printf("%lf ",x[h][i]);
}