三次样条插值算法

 1 #include<iostream>
 2 #include<math.h>
 3 using namespace std;
 4 float Hermite(int n, float xx)
 5 {
 6     float *x = new float[n+1],*ar = new float[n+1];
 7     float *y = new float[n+1],*br = new float[n+1];
 8     float *h = new float[n],m0,m1;
 9     float *a = new float[n+1], *b = new float[n+1];
10     float *m = new float[n+2],fx;
11     int i,choice=0;
12     for(i=0; i<n+1; i++)//输入xi，yi
13     {
14         cout<<"input x"<<i<<":";
15         cin>>x[i];
16         cout<<"input y"<<i<<":";
17         cin>>y[i];
18     }
19     for(i=0; i<n; i++)//计算hi
20     {
21         h[i] = x[i+1] - x[i];
22     }
23     for(i=1; i<n; i++)//i = 1,2,...n 计算α，β
24     {
25         ar[i] = h[i-1] / (h[i-1] + h[i]);
26         br[i] = 3 * ((1-ar[i])/h[i-1]*(y[i]-y[i-1])+ar[i]/h[i]*(y[i+1]-y[i]));
27     }
28     CHOICE:
29     cout<<"输入边界条件：\n1.第一种边界条件.\n2.第二种边界条件.\n"<<endl;
30     cin>>choice;//选择边界条件
31     switch (choice)
32     {
33     case 1:
34         cout<<"input M0:\n";
35         cin>>m0;
36         cout<<"input Mn:\n";
37         cin>>m1;
38         ar[0] = 0;        //第一种边界条件
39         br[0] = 2 * m0;
40         ar[n] = 1;
41         br[n] = 2 * m1;
42         break;
43     case 2:
44         ar[0] = 1;       //第二种边界条件
45         br[0] = 3 / h[0] * (y[1] - y[0]);
46         ar[n] = 0;
47         br[n] = 3 / h[n-1] * (y[n] - y[n-1]);
48         break;
49     default:
50         cout<<"input error!"<<endl;
51         goto CHOICE;
52         break;
53     }
54     //计算a,b
55     a[0] = -ar[0]/2;
56     b[0] = br[0]/2;
57     for(i=1; i<n+1; i++)
58     {
59         a[i] = -ar[i]/(2+(1-ar[i])*a[i-1]);
60         b[i] = (br[i]-(1-ar[i])*b[i-1])/(2+(1-ar[i])*a[i-1]);
61     }
62     //计算mi
63     m[n+1] = 0.0;
64     for(i=n; i>-1; i--)
65     {
66         m[i] = a[i] * m[i+1] + b[i];
67     }
68     if (xx < x[0] || xx > x[n])
69     {
70         cout<<"输入的数不在x取值范围内!"<<endl;
71         return 0.0;
72     }
73       for(i=0; i<n+1; i++)//输出mi
74     {
75         cout<<"m"<<i<<"="<<m[i]<<" ";
76         cout<<endl;
77     }
78     for(i=0; i<n; i++)
79     {
80         if(xx >= x[i] && xx < x[i+1])break;
81     }
82     //计算结果
83     fx = (1+2*(xx-x[i])/(x[i+1]-x[i]))*pow(((xx-x[i+1])/(x[i]-x[i+1])),2)*y[i] +
84          (1+2*(xx-x[i+1])/(x[i]-x[i+1]))*pow(((xx-x[i])/(x[i+1]-x[i])),2)*y[i+1] +
85          (xx-x[i])*pow(((xx-x[i+1])/(x[i]-x[i+1])),2)*m[i] +
86          (xx-x[i+1])*pow(((xx-x[i])/(x[i+1]-x[i])),2)*m[i+1];
87     cout<<"The answer is:\n"<<"s("<<xx<<")="<<fx<<endl;    //输出结果
88     return fx;
89 }
90 int main()
91 {
92     int n;
93     float xx;
94     cout<<"输入点的个数:\n";
95     cin>>n;
96     cout<<"输入x的值:\n";
97     cin>>xx;
98     Hermite(n-1,xx);
99 }