Akima 插值
附带的图片为运行结果
#include "stdio.h"
#include "math.h"
#include "interpolation.h"
void interpolation_akima(AKINTEP ap) {
int num,k,kk,m,l;
double pio,*mtr,*x,*y,u[5],p,q;
num=ap->n; k=ap->k;
pio=ap->t; mtr=ap->s;
x=ap->x; y=ap->y;
if (num<1) {
return;
}
else if (num==1) {
mtr[0]=mtr[4]=y[0];
return;
}
else if (num==2) {
mtr[0]=y[0];
mtr[1]=(y[1]-y[0])/(x[1]-x[0]);
if (k<0)
mtr[4]=(y[0]*(pio-x[1])-y[1]*(pio-x[0]))/(x[0]-x[1]);
return;
}
if (k<0) {
if (pio<=x[1]) kk=0;
else if (pio>=x[num-1]) kk=num-2;
else {
kk=1; m=num;
while (((kk-m)!=1)&&((kk-m)!=-1)) {
l=(kk+m)/2;
if (pio<x[l-1]) m=l;
else kk=l;
}
kk--;
}
}
else kk=k;
if (kk>=num-1) kk=num-2;
u[2]=(y[kk+1]-y[kk])/(x[kk+1]-x[kk]);
if (num==3) {
if (kk==0) {
u[3]=(y[2]-y[1])/(x[2]-x[1]);
u[4]=2.0*u[3]-u[2];
u[1]=2.0*u[2]-u[3];
u[0]=2.0*u[1]-u[2];
}
else {
u[1]=(y[1]-y[0])/(x[1]-x[0]);
u[0]=2.0*u[1]-u[2];
u[3]=2.0*u[2]-u[1];
u[4]=2.0*u[3]-u[2];
}
}
else {
if (kk<=1){
u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
if (kk==1) {
u[1]=(y[1]-y[0])/(x[1]-x[0]);
u[0]=2.0*u[1]-u[2];
if (num==4) u[4]=2.0*u[3]-u[2];
else u[4]=(y[4]-y[3])/(x[4]-x[3]);
}
else {
u[1]=2.0*u[2]-u[3];
u[0]=2.0*u[1]-u[2];
u[4]=(y[3]-y[2])/(x[3]-x[2]);
}
}
else if (kk>=(num-3)) {
u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
if (kk==(num-3)) {
u[3]=(y[num-1]-y[num-2])/(x[num-1]-x[num-2]);
u[4]=2.0*u[3]-u[2];
if (num==4) u[0]=2.0*u[1]-u[2];
else u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
}
else {
u[3]=2.0*u[2]-u[1];
u[4]=2.0*u[3]-u[2];
u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
}
}
else {
u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
u[4]=(y[kk+3]-y[kk+2])/(x[kk+3]-x[kk+2]);
}
}
mtr[0]=fabs(u[3]-u[2]);
mtr[1]=fabs(u[0]-u[1]);
if ((fabs(mtr[0])<0.0000001)&&(fabs(mtr[1])<0.0000001))
p=(u[1]+u[2])/2.0;
else p=(mtr[0]*u[1]+mtr[1]*u[2])/(mtr[0]+mtr[1]);
mtr[0]=fabs(u[3]-u[4]);
mtr[1]=fabs(u[2]-u[1]);
if ((fabs(mtr[0])<0.0000001)&&(fabs(mtr[1])<0.0000001))
q=(u[2]+u[3])/2.0;
else q=(mtr[0]*u[2]+mtr[1]*u[3])/(mtr[0]+mtr[1]);
mtr[0]=y[kk];
mtr[1]=p;
mtr[3]=x[kk+1]-x[kk];
mtr[2]=(3.0*u[2]-2.0*p-q)/mtr[3];
mtr[3]=(q+p-2.0*u[2])/(mtr[3]*mtr[3]);
if (k<0) {
p=pio-x[kk];
mtr[4]=mtr[0]+mtr[1]*p+mtr[2]*p*p+mtr[3]*p*p*p;
}
return;
}
main()
{
double x[11]={3.0,5.0,8.0,13.0,17.0,25.0,27.0,29.0,31.0,35.0,39.0};
double y[11]={7.0,10.0,11.0,17.0,23.0,18.0,13.0,6.0,3.0,1.0,0.0};
AKINTE aa= {11, x, y, -1, 14.0, {0}};
AKINTE ab= {11, x, y, -1, 28.0, {0}};
printf("\n");
interpolation_akima(&aa);
printf("x=%6.3f, f(x)=%e\n",aa.t, aa.s[4]);
printf("mtr0=%e, mtr1=%e, mtr2=%e, mtr3=%e\n",aa.s[0],aa.s[1],aa.s[2],aa.s[3]);
printf("\n");
interpolation_akima(&ab);
printf("x=%6.3f, f(x)=%e\n",ab.t, ab.s[4]);
printf("mtr0=%e, mtr1=%e, mtr2=%e, mtr3=%e\n",ab.s[0],ab.s[1],ab.s[2],ab.s[3]);
printf("\n");
}
三次样条插值的实现
1、程序比较简单的：
#include<iostream>
#include<iomanip>
using namespace std;
const int MAX = 50;
float x[MAX], y[MAX], h[MAX];
float c[MAX], a[MAX], fxym[MAX];
float f(int x1, int x2, int x3){
float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);
float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
return (a - b)/(x[x3] - x[x1]);
} //求差分
void cal_m(int n){ //用追赶法求解出弯矩向量M……
float B[MAX];
B[0] = c[0] / 2;
for(int i = 1; i < n; i++)
B[i] = c[i] / (2 - a[i]*B[i-1]);
fxym[0] = fxym[0] / 2;
for(i = 1; i <= n; i++)
fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);
for(i = n-1; i >= 0; i--)
fxym[i] = fxym[i] - B[i]*fxym[i+1];
}
void printout(int n);
int main(){
int n,i; char ch;
do{
cout<<"Please put in the number of the dots:";
cin>>n;
for(i = 0; i <= n; i++){
cout<<"Please put in X"<<i<<':';
cin>>x[i]; //cout<<endl;
cout<<"Please put in Y"<<i<<':';
cin>>y[i]; //cout<<endl;
}
for(i = 0; i < n; i++) //求 步长
h[i] = x[i+1] - x[i];
cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";
int t;
float f0, f1;
cin>>t;
switch(t){
case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";
cin>>f0>>f1;
c[0] = 1; a[n] = 1;
fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];
break;
case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";
cin>>f0>>f1;
c[0] = a[n] = 0;
fxym[0] = 2*f0; fxym[n] = 2*f1;
break;
default:cout<<"不可用\n";//待定
};//switch
for(i = 1; i < n; i++)
fxym[i] = 6 * f(i-1, i, i+1);
for(i = 1; i < n; i++){
a[i] = h[i-1] / (h[i] + h[i-1]);
c[i] = 1 - a[i];
}
a[n] = h[n-1] / (h[n-1] + h[n]);
cal_m(n);
cout<<"\n输出三次样条插值函数：\n";
printout(n);
cout<<"Do you to have anther try ? y/n :";
cin>>ch;
}while(ch == 'y' || ch == 'Y');
return 0;
}
void printout(int n){
cout<<setprecision(6);
for(int i = 0; i < n; i++){
cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";
/*
cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "
<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "
<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";
cout<<endl;*/
float t = fxym[i]/(6*h[i]);
if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";
else cout<<-t<<"*(x - "<<x[i+1]<<")^3";
t = fxym[i+1]/(6*h[i]);
if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";
else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";
cout<<"\n\t";
t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];
if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";
else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";
t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];
if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";
else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";
cout<<endl<<endl;
}
cout<<endl;
}
2、程序比较复杂的：
（程序前面的01.，02.，03.等等为语句编号，实际应用时请一一删除）01./*=======================================================================*/
02.#include <stdio.h>
03.////////////////////////////////////////////////////////////////////////////////
04.#define MAXNUM 50 //定义样条数据区间个数最多为50个
05.typedef struct SPLINE //定义样条结构体，用于存储一条样条的所有信息
06.{ //初始化数据输入
07. float x[MAXNUM+1]; //存储样条上的点的x坐标，最多51个点
08. float y[MAXNUM+1]; //存储样条上的点的y坐标，最多51个点
09. unsigned int point_num; //存储样条上的实际的 点 的个数
10. float begin_k1; //开始点的一阶导数信息
11. float end_k1; //终止点的一阶导数信息
12. //float begin_k2; //开始点的二阶导数信息
13. //float end_k2; //终止点的二阶导数信息
14. //计算所得的样条函数S(x)
15. float k1[MAXNUM+1]; //所有点的一阶导数信息
16. float k2[MAXNUM+1]; //所有点的二阶导数信息
17. //51个点之间有50个段，func[]存储每段的函数系数
18. float a3[MAXNUM],a1[MAXNUM];
19. float b3[MAXNUM],b1[MAXNUM];
20. //分段函数的形式为 Si(x) = a3[i] * {x(i+1) - x}^3 + a1[i] * {x(i+1) - x} +
21. // b3[i] * {x - x(i)}^3 + b1[i] * {x - x(i)}
22. //xi为x[i]的值，xi_1为x[i+1]的值
23.}SPLINE,*pSPLINE;
24.typedef int RESULT; //返回函数执行的结果状态，下面为具体的返回选项
25.#ifndef TRUE
26.#define TRUE 1
27.#endif
28.#ifndef FALSE
29.#define FALSE -1
30.#endif
31.#ifndef NULL
32.#define NULL 0
33.#endif
34.#ifndef ERR
35.#define ERR -2
36.#endif
37.//////////////////////////////////////////////////////////////////////////////////
38./*===============================================================================
39.*** 函数名称： Spline3()
40.*** 功能说明： 完成三次样条差值，其中使用追赶法求解M矩阵
41.*** 入口参数： (pSPLINE)pLine 样条结构体指针pLine中的x[],y[],num,begin_k1,end_k1
42.*** 出口参数： (pSPLINE)pLine 样条结构体指针pLine中的函数参数
43.*** 返回参数： 返回程序执行结果的状态TRUE or FALSE
44.================================================================================*/
45.RESULT Spline3(pSPLINE pLine)
46.{
47. float H[MAXNUM] = {0}; //小区间的步长
48. float Fi[MAXNUM] = {0}; //中间量
49. float U[MAXNUM+1] = {0}; //中间量
50. float A[MAXNUM+1] = {0}; //中间量
51. float D[MAXNUM+1] = {0}; //中间量
52. float M[MAXNUM+1] = {0}; //M矩阵
53. float B[MAXNUM+1] = {0}; //追赶法中间量
54. float Y[MAXNUM+1] = {0}; //追赶法中间变量
55. int i = 0;
56. ////////////////////////////////////////计算中间参数
57. if((pLine->point_num < 3) || (pLine->point_num > MAXNUM + 1))
58. {
59. return ERR; //输入数据点个数太少或太多
60. }
61. for(i = 0;i <= pLine->point_num - 2;i++)
62. { //求H[i]
63. H[i] = pLine->x[i+1] - pLine->x[i];
64. Fi[i] = (pLine->y[i+1] - pLine->y[i]) / H[i]; //求F[x(i),x(i+1)]
65. }
66. for(i = 1;i <= pLine->point_num - 2;i++)
67. { //求U[i]和A[i]和D[i]
68. U[i] = H[i-1] / (H[i-1] + H[i]);
69. A[i] = H[i] / (H[i-1] + H[i]);
70. D[i] = 6 * (Fi[i] - Fi[i-1]) / (H[i-1] + H[i]);
71. }
72. //若边界条件为1号条件，则
73. U[i] = 1;
74. A[0] = 1;
75. D[0] = 6 * (Fi[0] - pLine->begin_k1) / H[0];
76. D[i] = 6 * (pLine->end_k1 - Fi[i-1]) / H[i-1];
77. //若边界条件为2号条件，则
78. //U[i] = 0;
79. //A[0] = 0;
80. //D[0] = 2 * begin_k2;
81. //D[i] = 2 * end_k2;
82. /////////////////////////////////////////追赶法求解M矩阵
83. B[0] = A[0] / 2;
84. for(i = 1;i <= pLine->point_num - 2;i++)
85. {
86. B[i] = A[i] / (2 - U[i] * B[i-1]);
87. }
88. Y[0] = D[0] / 2;
89. for(i = 1;i <= pLine->point_num - 1;i++)
90. {
91. Y[i] = (D[i] - U[i] * Y[i-1]) / (2 - U[i] * B[i-1]);
92. }
93. M[pLine->point_num - 1] = Y[pLine->point_num - 1];
94. for(i = pLine->point_num - 1;i > 0;i--)
95. {
96. M[i-1] = Y[i-1] - B[i-1] * M[i];
97. }
98. //////////////////////////////////////////计算方程组最终结果
99. for(i = 0;i <= pLine->point_num - 2;i++)
100. {
101. pLine->a3[i] = M[i] / (6 * H[i]);
102. pLine->a1[i] = (pLine->y[i] - M[i] * H[i] * H[i] / 6) / H[i];
103. pLine->b3[i] = M[i+1] / (6 * H[i]);
104. pLine->b1[i] = (pLine->y[i+1] - M[i+1] * H[i] * H[i] / 6) /H[i];
105. }
106. return TRUE;
107.}
108.//////////////////////////////////////////////////////////////////////////////////
109.SPLINE line1;
110.pSPLINE pLine1 = &line1;
111.//////////////////////////////////////////////////////////////////////////////////
112.main()
113.{
114. line1.x[0] = 27.7;
115. line1.x[1] = 28;
116. line1.x[2] = 29;
117. line1.x[3] = 30;
118. line1.y[0] = 4.1;
119. line1.y[1] = 4.3;
120. line1.y[2] = 4.1;
121. line1.y[3] = 3.0;
122. line1.point_num = 4;
123. line1.begin_k1 = 3.0;
124. line1.end_k1 = -4.0;
125. Spline3(pLine1);
126. return 0;
127.}
128.//////////////////////////////////////////////////////////////////////////////////
