大學時的幾道數值演算法實現(code)

  1 多項式ax2+bx+c求解

/*        author: denny wqf363@hotmail.com       date：2004-12-19 16:04       e.g: x2+x+1       input: 1 1 1       output:                  two conjugated complex number roots:              x1=-0.50+ 0.87i              x2=-0.50- 0.87i       arithmetic：迭代求解*/#include <stdio.h>#include <math.h>       /*為使用開方函數sqrt( )包含數學函數庫*/int main(){       float a,b,c,d,e,f,x1,x2;       printf("input equation’s 3 quotieties :");       scanf("%f%f%f",&a,&b,&c);       if (fabs(a)<1e-4)         /*注意實數的精確比較方法*/              if (fabs(b)<1e-4)                     printf(" equation is invalid");              else                      printf("sigle root: x=%5.2f",-c/b);       else        {               d=-b/(2*a);              e=b*b-4*a*c;              if (e>0)                        {                      x1=d+sqrt(e)/(2*a);                     x2=d-sqrt(e)/(2*a);                     printf("two unequal real roots:/n x1=%5.2f/n x2=%5.2f/n",x1,x2);              }              else   if (fabs(e)<1e-4)                     printf("two equal roots : x=%5.2f/n",d);              else                 {                         f=sqrt(-e)/(2*a);                     printf("two conjugated complex number roots:/n");                     printf("x1=%5.2f+%5.2fi/n",d,f);                     printf("x2=%5.2f-%5.2fi/n",d,f);              }       }       return 0;} 2 迭代法求方程的根/*       Source of subject: 計算方法 page200 例3       author: denny wqf363@hotmail.com       date: 2004-12-19 16:04       E.g: x=log10(x+2)            x0=1       Output:0.3758       Algorithm: 使用迭代法求方程的根　      x_k+1=z(x_k) k=0,1...*/#include <math.h>#include <stdio.h>double z(double x){       return (log10(x+2));}iteration(){       double x,temp;       int k=0;       x=temp=1;       do       {              temp=x; /*save xk for compare*/              x=z(x);              k++;              printf("the %d time is: %f/n",k,x);       }while(fabs(x-temp)>=0.0001);}main(){       iteration();       getchar();    /*return until enter key!*/       return 0;} 3 經典四級四階R_K格式/*        source：1.3 page21 例1.3   經典四級四階R_K格式       author: denny wqf363@hotmail.com       date：2004-12-19        e.g: y'=y-2x/y,x∈[0,0.6]            y(0)=1       input: please enter the x0,y0: 1 -2       output:                  The table (x,y) is:               (1.200000,0.381267)               (1.400000,0.729676)               (1.600000,1.051183)               (1.800000,1.350665)               (2.000000,1.632115)               (2.200000,1.898800)               (2.400000,2.153398)               (2.600000,2.398099)               (2.800000,2.634697)               (3.000000,2.864661)       arithmetic：使用經典四階龍格-庫塔演算法進行高精度求解　      y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6       K1=h*f(xi,yi)      K2=h*f(xi+h/2,yi+K1/2)       K3=h*f(xi+h/2,yi+K2/2)      K4=h*f(xi+h,yi+K3)*/#include"stdio.h"#include"math.h"float f(float x,float y){       return (-y+x+1);}R_K44(){       float x0,y0,x,y;       float k1,k2,k3,k4;       float h=0.2;       printf("please enter the x0,y0: ");       scanf("%f,%f",&x0,&y0);       printf("The table (x,y) is:/n");       y=y0;                   /*賦初值*/       for(x=x0;x<3;x+=h)    {       k1=f(x,y);    k2=f(x+h/2,y+h*k1/2);    k3=f(x+h/2,y+h*k2/2);    k4=f(x+h,y+h*k3);    y+=h/6*(k1+2*k2+2*k3+k4);           /*y_x=0.5*((x+h)*(x+h)+2)*exp(-(x+h));*/ /*精確值*/           printf(" (%f,%f)/n ",x+h,y);     }}main(){       R_K44();}  4 Runge_Kutta_Fehlberg法#include <stdio.h>#include <math.h>/* page25 例1.4 Runge_Kutta_Fehlberg法*/main(){ /*輸入*/ double f(double,double); double eps,q,r,k1,k2,k3,k4,k5,k6; double x,x0,xn;          /*定義X的起點，終點*/ double y,y0,y_x,y_n,y__n;     /*定義Y的初值，精確解，Y_N是r_k四級四階解，y__n是r_K_f解*/ double h,hmax,hmin;       /*定義最大步長，最小步長*/ scanf("%f%f%f%f%f%f",&x0,&xn,&y0,&eps,&hmax,&hmin); x=x0;y=y0;h=hmax; printf("初值：(%f,%f)",x,y);  for(x=x0;x<xn;x+=h) {   k1=h*f(x,y);   k2=h*f(x+h/4,y+k1/4);   k3=h*f(x+3/8*h,y+3/32*k1+9/32*k2);   k4=h*f(x+12/13*h,y+1932/2197*k1-7200/2197*k2+7296/2197*k3);   k5=h*f(x+h,y+439/216*k1-8*k2+3680/513*k3-845/4104*k4);   k6=h*f(x+h/2,y-8/27*k1+2*k2-3544/2565*k3+1859/4104*k4-11/40*k5);   y_n+=25/216*k1+1408/2565*k3+2197/4104*k4-1/5*k5;   y__n+=16/135*k1+6656/12825*k3+28561/56430*k4-9/50*k5+2/55*k6;   r=(y__n-y_n)/h;          /*誤差R*/   q=0.84*pow(eps/r,0.25);      /*控制步長Q*/    if(r<=eps)    /*滿足精度要求，用r_k四級四階公式*/   {printf("(%f,%f,%f)",x,y_n,h);   /*確定最佳步長*/    if(q<=0.1)h=0.1*h;      else if(q>=4.0)h=4*h;        else h=q*h;    if(h>hmax)h=hmax;    if(h<hmin)printf("fail");   }   y_x=(x+h)+exp(-(x+h));   printf(" (%1.8f,%1.8f,%1.8f,%1.8f) ",x+h,y_x,y__n,y__n-y_x); /*x,精確解，R-K_F解，誤差*/ }}double f(double x,double y){double z; z=-y+x+1; return (z);} 5 主元高斯消去法#include <stdio.h>#include <math.h>#define N 50main(){  float b[N*(N+1)],a[N][N+1],max,temp,q; int p[N],i,j,n,k,row,col,l=0; printf("輸入方程組的未知數個數:/n"); scanf("%d",&n); printf("輸入方程組的增廣矩陣:/n"); for(i=0;i<n*(n+1);i++) scanf("%f",&b[i]); for(i=0;i<n;i++)   for(j=0;j<=n;j++)    a[i][j]=b[i*(n+1)+j]; for(i=0;i<n;i++)   p[i]=0; for(k=0;k<n;k++) { max=0;    for(i=0;i<n;i++)    { if(p[i]!=1)      {for(j=0;j<n;j++)       { if(p[j]!=1)        {if(fabs(max)<fabs(a[i][j]))        { max=a[i][j];row=i;col=j;}        }       }      }    }    if(max==0){printf("No only one");l=1;}    if(row!=col)    for(i=0;i<=n;i++)    { temp=a[col][i];a[col][i]=a[row][i];      a[row][i]=temp;    }    q=1.0/a[col][col];    for(i=0;i<n+1;i++)     a[col][i]*=q;    for(i=0;i<n;i++)    { if(i!=col)      { temp=a[i][col];          for(j=0;j<n+1;j++)        a[i][j]-=(a[col][j]*temp);      }    }   p[col]=1; } if(l==0) printf("方程解為:/n"); for(i=0;i<n;i++)   printf("x%d=%f/n",i,a[i][n]);} /*運行測試結果:(方程1)輸入方程組的未知數個數:3輸入方程組的增廣矩陣:2 1 -3 54 5 1 31 2 1 0方程解為:x0=1.000000x1=0.000000x2=-1.000000 (方程2)輸入方程組的未知數個數:3輸入方程組的增廣矩陣:0.0001 0.5402 0.3425 0.88281.235 2.567 0.9750 4.7771.024 2.001 4.555 7.580方程解為:x0=1.000000x1=1.000000x2=1.000000*/