llf.cpp

运用C++编写的一个二次插值运算程序

#include "stdio.h"
#include "math.h"
#include "conio.h"
///////////////////////////////////////////////////////////////////////////////////////////
void main()
{
  double *area(double a1,double p,double a[3]);
  double f(double x);
  double ar,fr;
  double a1=10.0,p=0.01,e=0.00001;
  double pa[3];
  area(a1,p,pa);
  a1=pa[0];
  double a2=pa[1];
  double a3=pa[2];
  double f1=f(a1);
  double f2=f(a2);
  double f3=f(a3);
 do
   {
      ar=(a3*a3-a2*a2)*f1+(a1*a1-a3*a3)*f2+(a2*a2-a1*a1)*f3;
      ar=ar/2/((a3-a2)*f1+(a1-a3)*f2+(a2-a1)*f3);
      fr=f(ar);
      if(ar>a2)
          {
             if(fr>f2)
               { 
                  a3=ar; f3=fr;
               }
             else if(fr<f2)
               { 
                   a1=a2;  f1=f2;
                   a2=ar;  f2=fr;    
                }
             else
              { 
                a3=ar; a1=a2; a2=(a1+a3)/2;
                f1=f2; f3=fr; f2=f(a2);   
              }
          }
      else if(ar<a2)
          {
             if(fr>f2)
               { 
                  a1=ar; f1=fr;
               }
             else if(fr<f2)
               { 
                   a3=a2;  f3=f2;
                   a2=ar;  f2=fr;    
                }
             else
              { 
                a1=ar; a3=a2; a2=(a1+a3)/2;
                f1=f2; f3=fr; f2=f(a2);   
              }
          }
      if(fabs(a1-a3)<=e)  break;
   }while(1);
  if(f2<fr)
     { ar=a2; fr=f2;}
  printf("\n自变量x*=%f\n因变量f(x*)=%f",ar,fr);
}


/////////////////////////////////////////////////////////////////////////////////////////


double *area(double a1,double p,double a[3])
{
  double f(double x);
  double a2,f2,a3,f3,temp;
  double acc=0.00001;
  double f1=f(a1);
  //double storep=p;
  while(1)
     { 
        a2=a1+p;    f2=f(a2);
        if(f2>=f1)
          {
             if(fabs(f2-f1)<acc)
                p=p/2;
             else
                p=-p;
          } 
        else break;
     }
    
  while(1)
     { 
        a3=a2+p;    f3=f(a3);
        if(f2<=f3)  break;
        p=p*2;
        a1=a2;      f1=f2;
        a2=a3;      f2=f3;
     }

  if(a1>a3)
    { 
       temp=a1;  a1=a3;  a3=temp;
    }
  a[0]=a1; a[1]=a2;  a[2]=a3;
  return a;
}

/////////////////////////////////////////////////////////////////////////////////////////////

double f(double x)
{
   double y=/*x*x+2*x+3+100;//*/pow(fabs(x-1),1.5)+pow(fabs(x-1),2.7);
   return y;
}