三次样条函数.c

实现 了 数值分析 课程 中得 几个 算法

/* 这份源代码文件已被未注册的SourceFormatX格式化过 */
/* 如果您想不再添加此类信息，请您注册这个共享软件  */
/* 更多相关信息请访问网站: http://cn.textrush.com  */


///////////////////////////////////////////////////////////////////
//            Purpose: 给定M0,Mn值的三次样条插值多项式           //
///////////////////////////////////////////////////////////////////


/**
 *The rusult:
 *Input n value: 3
 *Now input the (x_i,y_i),i=0,...,3:
 *1.1 0.4 1.2 0.8 1.4 1.65 1.5 1.8
 *Now input the M[0] value: 0
 *Now input the M[n] value: 0
 *Now input the x value: 1.25
 *S(1.250000)=1.033594
 *请按任意键继续. . .
 */
#include <stdio.h>
#define MAX_N 20            //定义(x_i,y_i)的最大维数

typedef struct tagPOINT  //点的结构
{
  double x;
  double y;
} POINT;

int main()
{
  int n;
  int i, k;
  POINT points[MAX_N + 1];
  double h[MAX_N + 1], b[MAX_N + 1], c[MAX_N + 1], d[MAX_N + 1], M[MAX_N + 1];
  double u[MAX_N + 1], v[MAX_N + 1], y[MAX_N + 1];
  double x, p, q, S;
  printf("\nInput n value: ");
  scanf("%d", &n);
  if (n > MAX_N)
  {
    printf("The input n is larger than MAX_N, please redefine the MAX_N.\n");
    return 1;
  }
  if (n <= 0)
  {
    printf("Please input a number between 1 and %d.\n", MAX_N);
    return 1;
  }
  //输入被插值点(x_i,y_i),M0值和Mn值
  printf("Now input the (x_i,y_i),i=0,...,%d:\n", n);
  for (i = 0; i <= n; i++)
    scanf("%lf%lf", &points[i].x, &points[i].y);
  printf("Now input the M[0] value: ");
  scanf("%lf", &M[0]);
  printf("Now input the M[n] value: ");
  scanf("%lf", &M[n]);
  printf("Now input the x value: "); //输入计算三次样条插值函数的x值
  scanf("%lf", &x);
  if (x > points[n].x || x < points[0].x)
  {
    printf("Please input a number between %f and %f.\n", points[0].x,
      points[n].x);
    return 1;
  }
  //计算M关系式中各参数的值
  h[0] = points[1].x - points[0].x;
  for (i = 1; i < n; i++)
  {
    h[i] = points[i + 1].x - points[i].x;
    b[i] = h[i] / (h[i] + h[i - 1]);
    c[i] = 1-b[i];
    d[i] = 6 *((points[i + 1].y - points[i].y) / h[i] - (points[i].y - points[i
      - 1].y) / h[i - 1]) / (h[i] + h[i - 1]);
  }
  //用追赶法计算Mi,i=1,...,n-1
  d[1] -= c[1] *M[0];
  d[n - 1] -= b[n - 1] *M[n];
  b[n - 1] = 0;
  c[1] = 0;
  v[0] = 0;
  for (i = 1; i < n; i++)
  {
    u[i] = 2-c[i] *v[i - 1];
    v[i] = b[i] / u[i];
    y[i] = (d[i] - c[i] *y[i - 1]) / u[i];
  }
  for (i = 1; i < n; i++)
  {
    M[n - i] = y[n - i] - v[n - i] *M[n - i + 1];
  }
  //计算三次样条插值函数在x处的值
  k = 0;
  while (x >= points[k].x)
    k++;
  k = k - 1;
  p = points[k + 1].x - x;
  q = x - points[k].x;
  S = (p *p * p * M[k] + q * q * q * M[k + 1]) / (6 *h[k]) + (p *points[k].y +
    q * points[k + 1].y) / h[k] - h[k]*(p *M[k] + q * M[k + 1]) / 6;
  printf("S(%f)=%f\n", x, S);

  //输出
  system("pause");
}