insert.h.cpp

数值计算利用filter分解法来求解线性方程组,简明方便

//二维三点插值法
  double slq3(double *x,double *y,double *z,int n,int m,double u,double v)
  { int nn,mm,ip,iq,i,j,k,l;
    double b[3],h,w;
    nn=3;
    if (n<=3) { ip=0;  nn=n;}
    else if (u<=x[1]) ip=0;
    else if (u>=x[n-2]) ip=n-3;
    else
      { i=1; j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (u<x[l-1]) j=l;
            else i=l;
          }
        if (fabs(u-x[i-1])<fabs(u-x[j-1])) ip=i-2;
        else ip=i-1;
      }
    mm=3;
    if (m<=3) { iq=0; mm=m;}
    else if (v<=y[1]) iq=0;
    else if (v>=y[m-2]) iq=m-3;
    else
      { i=1; j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (v<y[l-1]) j=l;
            else i=l;
          }
        if (fabs(v-y[i-1])<fabs(v-y[j-1])) iq=i-2;
        else iq=i-1;
      }
    for (i=0;i<=nn-1;i++)
      { b[i]=0.0;
        for (j=0;j<=mm-1;j++)
          { k=m*(ip+i)+(iq+j);
            h=z[k];
            for (k=0;k<=mm-1;k++)
              if (k!=j)
                h=h*(v-y[iq+k])/(y[iq+j]-y[iq+k]);
            b[i]=b[i]+h;
          }
      }
    w=0.0;
    for (i=0;i<=nn-1;i++)
      { h=b[i];
        for (j=0;j<=nn-1;j++)
          if (j!=i)
            h=h*(u-x[ip+j])/(x[ip+i]-x[ip+j]);
        w=w+h;
      }
    return(w);
  }
  //////////////////////////////////
  //全区插值
double slgq(double *x,double *y,double *z,int n,int m,double u,double v)
  { int ip,ipp,i,j,l,iq,iqq,k;
    double h,w,b[10];
    if (u<=x[0]) { ip=1; ipp=4;}
    else if (u>=x[n-1]) { ip=n-3; ipp=n;}
    else
      { i=1; j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (u<x[l-1]) j=l;
            else i=l;
          }
        ip=i-3; ipp=i+4;
      }
    if (ip<1) ip=1;
    if (ipp>n) ipp=n;
    if (v<=y[0]) { iq=1; iqq=4;}
    else if (v>=y[m-1]) { iq=m-3; iqq=m;}
    else
      { i=1; j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (v<y[l-1]) j=l;
            else i=l;
          }
        iq=i-3; iqq=i+4;
      }
    if (iq<1) iq=1;
    if (iqq>m) iqq=m;
    for (i=ip-1;i<=ipp-1;i++)
      { b[i-ip+1]=0.0;
        for (j=iq-1;j<=iqq-1;j++)
          { h=z[m*i+j];
            for (k=iq-1;k<=iqq-1;k++)
              if (k!=j) h=h*(v-y[k])/(y[j]-y[k]);
            b[i-ip+1]=b[i-ip+1]+h;
          }
      }
    w=0.0;
    for (i=ip-1;i<=ipp-1;i++)
      { h=b[i-ip+1];
        for (j=ip-1;j<=ipp-1;j++)
          if (j!=i) h=h*(u-x[j])/(x[i]-x[j]);
        w=w+h;
      }
    return(w);
  }