利用三次样条插值进行数据拟合
0 引 言
三次样条插值以构造简单,使用方便,拟合准确,具有“保凸”的重要性质等特点成为了常用的插值方法。一般三次样条插值解算过程中通过追赶法求解三弯矩阵,但使用计算机求解时会表现出解的精度不高的问题,导致其计算结果无法应用到工程实践之中。因此需要找出一种提高解精度的方法。
1 基本概念
三次样条函数的定义:在区间内对于给定的函数值 ,其中,如果函数满足条件:
(1) 在每个子区间,上都是不高于三次的多项式;
(2) 、、在上都连续;
(3) , 。
则称为函数关于节点的三次样条函数。
想要求解三次样条插值函数,只需在每个子区间上确定一个三次多项式共有4个系数,确定它们需要 4n 个条件,因此要完全确定共需 4n 个条件。由所满足的条件 (1)、(2)、(3),可确定个条件,仍然缺少两个条件。这两个条件通常由实际问题对三次样条插值函数在端点的状态要求给出,也称之为边界条件,常见的边界条件有:
1)夹持边界条件(Clamped Spline):给定两端点的一阶导数值,即,;
2)自然边界条件(Natural Spline):使两端点的二阶导数值为零,即;
3)非扭结边界条件(Not-A-Knot Spline):强制第一个插值点的三阶导数值等于第二个点的三阶导数值,最后第一个点的三阶导数值等于最后第二个点的三阶导数值,即,。
2 计算方法
设三次样条函数,(0)
,,
,
由三次样条函数定义(1)(2)(3)可得:
, (1)
如下构造式(1)矩阵:
(2)
由式(1)可知:
,,,,
(3)
1)在夹持边界条件时,
,,,;
,,,;
2)在自然边界条件时,
, ,,;
,,,;
3)在非扭结边界条件时,
,,,;
,,,;
由n个未知数的非齐次方程组有惟一解的充分必要条件是,可知矩阵方程(2)在以上三种情况下都有惟一解。
对矩阵方程(2)采用高斯列主元消去法即可求解得出。
最后,代入式(0)可以得出:
,,,,
3 应用算例
有点集,在非扭结边界条件下进行插值。同时使用Matlab R2010a和文章所述方法进行插值计算,对比计算结果。
表1.Matlab R2010a计算结果
插值 方程 |
方程系数 |
|||
a | b | c | d | |
0.083 333 333 333 333 2 |
-0.5 |
1.416 666 666 666 67 |
1 |
|
0.083 333 333 333 333 3 |
-0.25 |
0.666 666 666 666 667 |
2 |
|
0.083 333 333 333 333 3 |
0.25 |
0.666 666 666 666 667 |
3 |
表2.文章所述方法计算结果
插值 方程 |
方程系数 |
|||
a | b | c | d | |
0.083 333 333 333 333 3 |
-0.5 |
1.416 666 666 666 665 |
1 |
|
0.083 333 333 333 333 3 |
-0.25 |
0.666 666 666 666 666 |
2 |
|
0.083 333 333 333 333 3 |
0.25 |
0.666 666 666 666 667 |
3 |
由表1、表2可知文章所述方法计算结果与Matlab R2010a计算一致,其他几种边界条件下的插值结论均一致,这里不在复述。
4 结 论
利用上述方法在传统的三弯矩阵中增加两个新元素,使得三次样条插值的求解过程可以在不增加计算复杂度的前提下提高计算结果精度。并且可以实现使用一种算法来计算三种常用边界条件下的样条插值函数。这样利于编制计算机程序来自动求解多种不同边界的三次样条函数。
1 /**********************************************************
2 ** 三次样条插值函数.C
3 ** 进行自然边界,夹持边界,非扭结边界条件下的插值
4 **
5 ** 2012-05-25 将函数从原来的需要至少4个插值点才能计算
6 ** 扩展成只需要2个插值点就可以完成计算
7 ** 其他一些改变
8 ** 2008-12-27 创建函数
9 **********************************************************/
10 #include <stdio.h>
11 #include <stdlib.h>
12
13 #define ABS(p) ( ((p) > 0) ? (p) : -(p) )
14 #define SWAP(x, y, temp) (temp) = (y); (y) = (x); (x) = (temp);
15
16 typedef enum _condition
17 {
18 NATURAL, // 自然边界
19 CLAMPED, // 夹持边界
20 NOTAKNOT // 非扭结边界
21 }condition;
22
23 typedef struct _coefficient
24 {
25 double a3;
26 double b2;
27 double c1;
28 double d0;
29 }coefficient;
30
31 typedef struct _point
32 {
33 coefficient *coe; // 插值结果系数矩阵
34 double *xCoordinate; // 插值点横坐标
35 double *yCoordinate; // 插值点纵坐标
36 double f0; // 在夹持条件下的最左边点的导数值
37 double fn; // 在夹持条件下的最右边点的导数值
38 int num; // 插值点数
39 condition con; // 边界条件
40 }point;
41
42
43 /*
44 ** 资源不足时函数返回 -1
45 ** 插值点数小于2时,函数返回 -2
46 ** 计算正确完成函数返回插值点的数量 n
47 */
48 int spline(point *point)
49 {
50 double *x = (*point).xCoordinate, *y = (*point).yCoordinate;
51 int n = (*point).num - 1; // 循环上限
52 int i; // 循环辅助变量
53 coefficient *coe = (*point).coe;
54 condition con = (*point).con;
55 double *h, *d;
56 double *a, *b, *c, *f, *m;
57 double temp;
58
59 if (n < 1)
60 {return -2;}
61
62 h = (double *)malloc( (6 * n + 4) * sizeof(double) ); /* 0, 1,...,(n-1) */
63 if (h == NULL)
64 {return -1;}
65 d = h + n; /* 0, 1,...,(n-1) */
66 a = d + n; /* 特别使用,1,...,(n-1), n */
67 b = a + (n + 1); /* 0,1,...,(n-1), n */
68 c = b + (n + 1); /* 0,1,...,(n-1),特别使用 */
69 f = c + (n + 1); /* 0,1,...,(n-1), n */
70 m = b;
71
72
73 /* 计算 h[] d[] */
74 for (i = 0; i < n; i++)
75 {
76 h[i] = x[i + 1] - x[i];
77 d[i] = (y[i + 1] - y[i]) / h[i];
78 /* printf("%f\t%f\n", h[i], d[i]); */
79 }
80 /**********************
81 ** 初始化系数增广矩阵
82 **********************/
83 a[0] = (*point).f0;
84 c[n] = (*point).fn;
85 /* 计算 a[] b[] c[] f[] */
86 switch(con)
87 {
88 case NATURAL:
89 f[0] = 0;
90 f[n] = 0;
91 a[0] = 0;
92 c[n] = 0;
93 c[0] = 0;
94 a[n] = 0;
95 b[0] = 1;
96 b[n] = 1;
97 break;
98
99 case CLAMPED:
100 f[0] = 6 * (d[0] - a[0]);
101 f[n] = 6 * (c[n] - d[n - 1]);
102 a[0] = 0;
103 c[n] = 0;
104 c[0] = h[0];
105 a[n] = h[n - 1];
106 b[0] = 2 * h[0];
107 b[n] = 2 * h[n - 1];
108 break;
109
110 case NOTAKNOT:
111 f[0] = 0;
112 f[n] = 0;
113 a[0] = h[0];
114 c[n] = h[n - 1];
115 c[0] = -(h[0] + h[1]);
116 a[n] = -(h[n - 2] + h[n - 1]);
117 b[0] = h[1];
118 b[n] = h[n - 2];
119 break;
120 }
121
122 for (i = 1; i < n; i++)
123 {
124 a[i] = h[i - 1];
125 b[i] = 2 * (h[i - 1] + h[i]);
126 c[i] = h[i];
127 f[i] = 6 * (d[i] - d[i - 1]);
128 }
129 /* for (i = 0; i < n+1; i++) printf("%f\n", c[i]); */
130
131 /*************
132 ** 高斯消元
133 *************/
134 if (n > 2)
135 {
136 /* 第0行到第(n-3)行的消元 */
137 for(i = 0; i <= n - 3; i++)
138 {
139 /* 选主元 */
140 if ( ABS(a[i + 1]) > ABS(b[i]) )
141 {
142 SWAP(a[i + 1], b[i], temp);
143 SWAP(b[i + 1], c[i], temp);
144 SWAP(c[i + 1], a[i], temp);
145 SWAP(f[i + 1], f[i], temp);
146 }
147 temp = a[i + 1] / b[i];
148 a[i + 1] = 0;
149 b[i + 1] = b[i + 1] - temp * c[i];
150 c[i + 1] = c[i + 1] - temp * a[i];
151 f[i + 1] = f[i + 1] - temp * f[i];
152 }
153 }
154 if (n >= 2)
155 {
156 /* 最后3行的消元 */
157 /* 第(n-2)行的消元 */
158 /* 选主元 */
159 if ( ABS(a[n - 1]) > ABS(b[n - 2]) )
160 {
161 SWAP(a[n - 1], b[n - 2], temp);
162 SWAP(b[n - 1], c[n - 2], temp);
163 SWAP(c[n - 1], a[n - 2], temp);
164 SWAP(f[n - 1], f[n - 2], temp);
165 }
166 /* 选主元 */
167 if ( ABS(c[n]) > ABS(b[n - 2]) )
168 {
169 SWAP(c[n], b[n - 2], temp);
170 SWAP(a[n], c[n - 2], temp);
171 SWAP(b[n], a[n - 2], temp);
172 SWAP(f[n], f[n - 2], temp);
173 }
174 /* 消元 */
175 temp = a[n - 1] / b[n - 2];
176 a[n - 1] = 0;
177 b[n - 1] = b[n - 1] - temp * c[n - 2];
178 c[n - 1] = c[n - 1] - temp * a[n - 2];
179 f[n - 1] = f[n - 1] - temp * f[n - 2];
180 /* 消元 */
181 temp = c[n] / b[n - 2];
182 c[n] = 0;
183 a[n] = a[n] - temp * c[n - 2];
184 b[n] = b[n] - temp * a[n - 2];
185 f[n] = f[n] - temp * f[n - 2];
186 }
187 /* 最后2行的消元 */
188 /* 第(n-1)行的消元 */
189 /* 选主元 */
190 if ( ABS(a[n]) > ABS(b[n - 1]) )
191 {
192 SWAP(a[n], b[n - 1], temp);
193 SWAP(b[n], c[n - 1], temp);
194 SWAP(f[n], f[n - 1], temp);
195 }
196 /* 消元 */
197 temp = a[n] / b[n-1];
198 a[n] = 0;
199 b[n] = b[n] - temp * c[n - 1];
200 f[n] = f[n] - temp * f[n - 1];
201
202 /******************
203 ** 回代求解 m[]
204 ******************/
205 m[n] = f[n] / b[n];
206 m[n - 1] = (f[n - 1] - c[n - 1] * m[n]) / b[n-1];
207 for (i = n - 2; i >= 0; i--)
208 {
209 m[i] = ( f[i] - (m[i + 2] * a[i] + m[i + 1] * c[i]) ) / b[i];
210 }
211 /* for (i = 0; i < n+1; i++) printf("%f\n", m[i]); */
212
213 /***********************
214 ** 计算插值函数所有参数
215 ***********************/
216 for (i = 0; i < n; i++)
217 {
218 coe[i].a3 = (m[i + 1] - m[i]) / (6 * h[i]);
219 coe[i].b2 = m[i] / 2;
220 coe[i].c1 = d[i] - (h[i] / 6) * (2 * m[i] + m[i + 1]);
221 coe[i].d0 = y[i];
222 }
223 free(h);
224 // 计算正确完成返回插值点的数量
225 return n + 1;
226 }
227
228
229 void main()
230 {
231 /* x[]内不能出现重复数据 */
232 double x[] = { 0, 0.5, 0.75, 1.25, 2.5, 5, 7.5, 10, 15,
233 20, 25, 30, 35, 40, 45, 50, 55, 60,
234 65, 70, 75, 80, 85, 90, 95, 100};
235 double y[] = { 0, 0.6107, 0.7424, 0.9470, 1.3074, 1.7773, 2.1,
236 2.3414, 2.6726, 2.8688, 2.9706, 3.0009, 2.9743, 2.9015,
237 2.7904, 2.6470, 2.4762, 2.2817, 2.0662, 1.8320, 1.5802,
238 1.3116, 1.0263, 0.7239, 0.4033, 0.0630};
239 int n = sizeof(x) / sizeof(double);
240 coefficient *coe;
241 int i;
242 point p;
243 coe = (coefficient *)malloc((n - 1) * sizeof(coefficient));
244 p.xCoordinate = x;
245 p.yCoordinate = y;
246 /* 下面两个参数f0和fn只在夹持条件下有效,其他条件下这两个值会被忽略 */
247 p.f0 = 0;// 在夹持条件下的左边点的导数值
248 p.fn = 0;// 在夹持条件下的右边点的导数值
249 p.num = n;
250 p.con = NOTAKNOT;
251 p.coe = coe;
252 spline(&p);
253
254 for (i = 0; i < n - 1; i++)
255 printf("%f\t%f\t%f\t%f\n", coe[i].a3, coe[i].b2, coe[i].c1, coe[i].d0);
256 free(coe);
257 }
[参 考 文 献]
[1] 李桂成.计算方法[M] 北京:电子工业出版社,2006.8
[2] 司守奎、孙玺菁.数学建模算法与应用[M] 北京:国防工业出版社,2011.8
[3] 同济大学数学系.线性代数[M] 北京:高等教育出版社,2007.5
[4] 于洋 袁健华 钱江 王美玲.新边界条件下的三次样条插值函数[J].软件 ,2016, (2)
[5] 何丽丽.三次样条插值--三转角方程的算法设计[J].湖南邮电职业技术学院学报 ,2014, (3)
[6] 俞海军 陈瑾怡.三种插值方法的研究与比较[J].河南科技 ,2013, (20)
[7] 乃明 胡兵 闵心畅.R-B方程样条有限元法[J].四川大学学报(自然科学版) ,2012, 49(5)
[8] 邢丽.三次样条插值端点约束条件的构造与Matlab实现[J].上海第二工业大学学报 ,2012, (4)
[9] 朱家明.两种三次样条插值建立方法及MATLAB运行时间比较[J].中国培训,2016, (10)
下一篇: C++中实现三次样条曲线
推荐阅读
-
利用曲线拟合进行插值和平滑处理
-
用三重矩方程和追赶法进行三次样条插值的数值方法实验(MATLAB 代码)
-
统计学习 04:假设检验(以 t 检验为例)和 P 值 - 要点 I. 假设检验的一般思路 假设检验 清楚你的问题是什么?期望得出什么结论? 例如,两种药物的疗效是否存在差异,自变量与因变量之间是否存在回归关系 .... 请始终牢记,假设检验回答的是是否存在某种关系的问题:它并不衡量这种关系有多大。 提出两种假设:零假设 (H0) 和备择假设 (H1) 零假设与备择假设相反,一般来说,研究的目的是证明原假设是错误的,即得出备择假设的结论。 例如,如果实验预期希望两种药物的疗效存在差异,那么 H0:μ1 - μ2 = 0;H1:μ1 - μ2 ≠ 0 H0:μ1-μ2 = 0 的一般形式称为双侧检验,而 >、<等零假设称为单侧检验。一般来说双侧检验更为常见,下面也主要介绍这种方法。 单尾或双尾测试 根据原始数据计算零假设概率分布的统计量(t 值、Z 值、F 值等)。 根据问题的性质选择合适的概率检验方法,从而计算出相应的统计量值;因此,不同情况的统计量值有不同的计算方法。 根据计算出的统计量值,利用统计软件,可以知道相应的 p 值是多少 也可以先确定一个合适的显著性水平(0.0.001....),并计算其临界值,再与我们计算出的统计量值进行比较,从而做出判断。 根据第四步的比较结果,如果 p 值小于预期的显著性水平(α,通常设定为 0.05),则认为该统计量远离原假设分布,属于小概率事件,则拒绝原假设,从而接受备择假设。 决定 要点 2:以 t 检验为例,演示上述假设检验思路。 t 检验基于 t 分布,常见的 t 检验有三种,如下图所示,但我认为第三种配对设计可能更常用(零假设:差异是否为零),下面介绍的例子就是一种配对设计 三次 t 检验 举例测量两组大鼠肝脏中维生素 A 的含量,比较两组大鼠维生素 A 含量是否有差异。数据如下 数据 (1) 预计两组大鼠的维生素 A 水平存在差异 (2) H0:μd=0,H1:μd≠0,α=0.05,双侧检验 (3) t 统计量的计算 配方 计算 上述程序计算的是*度为 7 的 t 分布情况下的 t 值。只需理解公式即可,不同的方法有不同的公式,这些交给统计软件即可。
-
Scipy:进行三次样条插值
-
MATLAB程序:用三次样条插值法进行数值分析
-
基于MATLAB的三维数据插值拟合与三次样条拟合算法(附完整代码)
-
使用MATLAB进行三次样条插值
-
如何使用Matlab进行三次样条插值
-
三次样条插值求拟合曲线 python 三次样条插值函数程序
-
利用三次样条插值进行数据拟合