计算方法中的一些程序

作者在 2008-05-02 12:10:17 发布以下内容

加速法

#include<stdio.h>
#include<math.h>
void main()
{
 float work_out(float x);
 float f1(float x);
 float a;
 printf("input origial:\n");
 scanf("%f",&a);
 f1(a);
 printf("%0.3f",work_out(a));
}
float f1(float x)
{
 
 x=(sin(x))/(x);
 return(x);
}

float work_out(float x)
{
 float x1,x2,x3;
  do{
  x1=f1(x);
  x2=f1(x1);
  x3=x2-((x2-x1)*(x2-x1))/(x2-2*x1+x);
 }while(fabs(x3-x2)>=1e-2);
  return(x3);
}

 

牛顿法

#include<stdio.h>
#include<math.h>
void main()
{
 float work_out(float x1);
 float x1,a;
 printf("输入x的一个值:\n");
 scanf("%f",&x1);
 a=work_out(x1);
 printf("输出计算结果:\n");
 printf("%0.4f",a);
}
float work_out(float x1)
{
 float x0,f0,f1;
 do {
  x0=x1;
  f0=x0-cos(x0);
  f1=1+sin(x0);
  x1=x0-(f0/f1);
 }while(fabs(x1-x0)>=1e-3);
 return(x1);
}

 

弦截法

#include<stdio.h>
#include<math.h>
float f(float x)
{
 float y;
 y=x*(x*(x-1))-1;
 return(y);
}
float xpoint(float x1,float x2)
{
 float y;
 y=(x1*f(x2)-x2*f(x1))/(f(x2)-f(x1));
 return(y);
}
float root(float x1,float x2)
{
 float x,y,y1;
 y1=f(x1);
 do{
  x=xpoint(x1,x2);
  y=f(x);
  if(y*y1>0)
  {
   y1=y;
   x1=x;
  }
  else
   x2=x;
 }while(fabs(y)>=1e-3);
  return(x);
}
void main()
{
 float x1,x2,f1,f2,x;
 do{
  printf("input x1 x2:\n");
  scanf("%f,%f",&x1,&x2);
  f1=f(x1);
  f2=f(x2);
 }
 while(f1*f2>=0);
 x=root(x1,x2);
 printf("A of equation is %0.4f\n",x);
}

 

迭代:
#include<stdio.h>
#include<conio.h>
#include<math.h>
#define n 3  
#define eps 0.5e-4  
#define Kmax 100
static double aa[n][n]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
 static double bb[n]={7.2,8.3,4.2};
void main()
{ int k,i,j;   double d,sum,norm;
  double a[n+1][n+1],b[n+1],x[n+1],y[n+1];
  for(i=1;i<=n;i++)
   { for(j=1;j<=n;j++)
   a[i][j]=aa[i-1][j-1];
     b[i]=bb[i-1];}
  for(i=1;i<=n;i++) x[i]=0;
  k=0;
  printf("\n  k=%2d    X=",k);
  for(i=1;i<=n;i++) printf(" %f  ",x[i]);
   do
    {  k++;  if(k>Kmax)
{ printf("\n The itetation failed !");break; }
       norm=0.0;
       for(i=1;i<=n;i++)
       { sum=0.0;
  for(j=1;j<=n;j++) if(j!=i) sum=sum+a[i][j]*x[j];
  y[i]=(b[i]-sum)/a[i][i];
  d=fabs(y[i]-x[i]);
  if(norm<d) norm=d;  
 }
      printf("\n  k=%2d    X=",k);
      for(i=1;i<=n;i++) printf(" %f  ",y[i]);
      for(i=1;i<=n;i++) x[i]=y[i];
     }
    while(norm>=eps); 
    if (norm<eps)
     { printf("\n\n The result is:\n");
      printf("\n  k=%d  ",k);
      for(i=1;i<=n;i++) printf("   x[%d]=%f",i,x[i]);
      }
   getch();
}

三角
#include<stdio.h>
#include<conio.h>
#include<math.h>
#define n 3
#define precision 1e-16
main()
{int i,j;double a[n+1][n+1],b[n+1];
 void LUdecomp();   void Solve();
 clrscr();
 printf("\n Please input the matrix A:\n");
 for(i=1;i<=n;i++)
  { printf("\n row %d:",i);
    for(j=1;j<=n;j++)  scanf("%lf,",&a[i][j]);
   }
 printf("\n Pleae input b[1]... b[n]:");
 for(i=1;i<=n;i++) scanf("%lf",&b[i]);
 LUdecomp(a);
 printf("\n LU decomposition of Doolittle :\n\n");
   for(i=1;i<=n;i++)
   { for(j=1;j<=n;j++) printf("%12f",a[i][j]);  printf("\n");
    }
 if(fabs(a[n][n])>precision)
  { Solve(a,b);  printf("\n");
    for(i=1;i<=n;i++)  printf("   x[%d]=%f",i,b[i]);
  }
  else  printf("\n LU decomposition method failed! \n");
getch();
  }
void LUdecomp(double a[][n+1])
 { int i,j,k,s;
   for(k=1;k<=n;k++)
    { for(j=k;j<=n;j++)  for(s=1;s<=k-1;s++)
 a[k][j]-=a[k][s]*a[s][j];
      if((fabs(a[k][k])<precision)&&(k<n))
     { printf("\n LU decomposition failed! \n"); exit(0);}
      for(i=k+1;i<=n;i++)
       { for(s=1;s<=k-1;s++)  a[i][k]-=a[i][s]*a[s][k];
      a[i][k]/=a[k][k];
 }
     }
  }
 void Solve(double a[][n+1],double b[])
  { int k,j;
    for(k=1;k<=n;k++) for(j=1;j<=k-1;j++)
   b[k]-=a[k][j]*b[j];
    for(k=n;k>=1;k--)
     { for(j=k+1;j<=n;j++) b[k]-=a[k][j]*b[j];
  b[k]/=a[k][k];
     }
  }

高斯
#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
#define n 3  
#define precision 1e-16  
static double aa[n][n+1]={{1,2,-1,3},{1,-1,5,0},{4,1,-2,2}};
main()
{ int i,j,det;  double a[n+1][n+2],x[n+1];
  int GaussElimination_ColumnSelect();
  clrscr();
  for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)
     a[i][j]=aa[i-1][j-1];
det=GaussElimination_ColumnSelect(a,x);
  if(det!=0)
    for(i=1;i<=n;i++)
     printf("\n  x[%d]=%f\n",i,x[i]); printf("\n");
  getch();
 }
int GaussElimination_ColumnSelect(double a[][n+2],double x[n+1])
 { int i,j,k,r;   double c;
   for(k=1;k<=n-1;k++)
    { r=k;
      for(i=k;i<=n;i++) 
 if(fabs(a[i][k])>fabs(a[r][k]))   r=i;
if(fabs(a[r][k])<precision)
{printf("\n det A = 0.Elimination  failed ! ");exit(0);}
      if(r!=k)
 {for(j=k;j<=n+1;j++) 
  { c=a[k][j];
a[k][j]=a[r][j]; a[r][j]=c;}
  }
      for(i=k+1;i<=n;i++)  
       { c=a[i][k]/a[k][k];
  for(j=k+1;j<=n+1;j++)
     a[i][j]=a[i][j]-c*a[k][j];
       }
     }
   if(fabs(a[n][n])<precision)
     {printf("\n det A = 0. Algorithm failed !");exit(0);}
   for(k=n;k>=1;k--) 
    { x[k]=a[k][n+1];
      for(j=k+1;j<=n;j++)
 x[k]=x[k]-a[k][j]*x[j];
      x[k]=x[k]/a[k][k];
      }
   return(1);
 }

 

 

文章评论,共0条
游客请输入验证码