22.cpp

三次样条插值计算方程组

// 22.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include<iostream.h>
#include<math.h>
#include<iomanip.h>

void splent (double xa[],double ya[],double y2a[],int n,
			 double &x,double&y)
 {
  int  klo,khi,k;
  double h ,a,b,aaa,bbb;
  klo=1;
  khi=n;
  loop:  if (khi-klo>1) 
        {
          k=(khi+klo)/2;
          if (xa[k]>x)
            khi=k;
          else
           {
        klo =k;
            }
      goto loop;
        }
        h=xa[khi]-xa[klo];
         if (h==0)
		 {
    cout<<"pause 'bad xa input'"<<endl;
    return;
  }
   a=(xa[khi]-x)/h;
   b=(x-xa[klo])/h;
   aaa=a*ya[klo]+b*ya[khi];
   bbb=(a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi];
   y=aaa+bbb*h*h/6.0;
}


void spline (double x[],double y[],int n ,double yp1,
			 double ypn,double y2[])
 {
   double u[100] ,aaa,sig,p,bbb,ccc,qn,un; 
   int  i,k;
   if (yp1>9.9e+29)
   {
   y2[1]=0;
   u[1]=0;
   }
  else
  {
   y2[1]=-0.5;
   aaa=(y[2]-y[1])/(x[2]-x[1]); 
   
   
   u[1]=(3/(x[2]-x[1]))*(aaa-yp1);
  }
  for (i=2;i<=n-1;i++)
  {
    sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p=sig*y2[i-1]+2.0;
    y2[i]=(sig-1.0)/p;
    aaa=(y[i+1]-y[i])/(x[i+1]-x[i]);
    bbb=(y[i]-y[i-1])/(x[i]-x[i-1]);
    ccc=x[i+1]-x[i-1];
    u[i]=(6.0*(aaa-bbb)/ccc-sig*u[i-1])/p;
  }
  if (ypn>9.9e+29)
  {
    qn=0.0;
    un=0.0;
  } 
  else
  {
     qn=0.5;
     aaa=ypn-(y[n]-y[n-1])/(x[n]-x[n-1]);
     un=(3.0/(x[n]-x[n-1]))*aaa;
  }
  y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  for (k=n-1;k>=1;k--) 
  y2[k]=y2[k]*y2[k+1]+u[k];
}


 void main( )
   {
   double x[21],y[21],y2[21],pi,yp1,ypn;
   int n,i;
   n=20;
   pi=3.1415926;
   //生成插值数组
   for (i=1;i<=20;i++)
  {
  x[i]=i*pi/n;
  y[i]=sin(x[i]);
  }
//用二次样条函数计算一阶导数
  yp1=cos(x[1]);
  ypn=cos(x[n]);
  cout<<endl;
  cout<<"second－dervative for sin(x) from 0 to pi "<<endl;
  spline(x,y,n,yp1,ypn,y2);
  //验证计算结果
  cout<<"         spline    actual"<<endl;
  cout<<"   angle  2nd derive    2nd derive" <<endl;
  for(i=1;i<=n;i++)
  {
  cout<<setprecision(6)<<setiosflags(ios::fixed);
  cout<<setw(14)<<x[i];
  cout<<setw(14)<<y2[i];
  cout<<setw(14)<<-sin(x[i])<<endl;
  }
}