三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。本篇介绍力求用容易理解的方式,介绍一下三次样条插值的原理,并附C语言的实现代码。

1. 三次样条曲线原理

假设有以下节点

image

 
 
1.1 定义

样条曲线image 是一个分段定义的公式。给定n+1个数据点,共有n个区间,三次样条方程满足以下条件:

a. 在每个分段区间image (i = 0, 1, …, n-1,x递增), image 都是一个三次多项式。

b. 满足image (i = 0, 1, …, n )

c. image ,导数image ,二阶导数image 在[a, b]区间都是连续的,即image曲线是光滑的。

所以n个三次多项式分段可以写作:

image ,i = 0, 1, …, n-1

其中ai, bi, ci, di代表4n个未知系数。

1.2 求解

已知:

a. n+1个数据点[xi, yi], i = 0, 1, …, n

b. 每一分段都是三次多项式函数曲线

c. 节点达到二阶连续

d. 左右两端点处特性(自然边界,固定边界,非节点边界)

根据定点,求出每段样条曲线方程中的系数,即可得到每段曲线的具体表达式。

 

插值和连续性:

image, 其中 i = 0, 1, …, n-1

微分连续性:

image , 其中 i = 0, 1, …, n-2

样条曲线的微分式:

imageimage

 

将步长 带入样条曲线的条件:

a. 由image (i = 0, 1, …, n-1)推出

image 

b. 由image (i = 0, 1, …, n-1)推出

image

c. 由 image (i &#

  • 19
    点赞
  • 70
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
#include <iostream> #include <cmath> #include <fstream> using namespace std; class scyt { float *x,*y,*d,*h,*u,*q,*a,*b,*c,*l,*r,*o,*M; int m; float y0,y3; public: scyt(); void qiudao(); void zgf(); void qiujie(); ~scyt(); }; void main() { scyt hello; hello.qiudao(); hello.zgf(); hello.qiujie(); } scyt::scyt() { ifstream fin("三次样条.txt"); for(float j;fin>>j;) { m=int(j); break; } x=new float[m]; y=new float[m]; d=new float[m]; h=new float[m-1]; u=new float[m-2]; q=new float[m-2]; a=new float[m-1]; b=new float[m]; c=new float[m-1]; l=new float[m]; r=new float[m-1];//此处的r为追赶法中的u; o=new float[m];//此处o为追赶法中的y M=new float[m];//此处M为追赶法中的x; int jishu=0; for(j;fin>>j;) { if(jishu<=m-1) x[jishu]=j; if(jishu>m-1&&jishu;<2*m) { y[jishu-m]=j; } if(jishu==2*m) { y0=j; } if(jishu==2*m+1) { y3=j; } jishu++; } fin.close(); } void scyt::qiudao() { for(int i=0;i<m-1;i++) { h[i]=x[i+1]-x[i]; } for(i=0;i<m-2;i++) { u[i]=h[i] / (h[i] + h[i+1]); } for(i=0;i<m-2;i++) { q[i]=1-u[i]; } d[0]=6/h[0]*((y[1]-y[0])/h[0]-y0); for(i=1;i<m-1;i++) { d[i]=6/(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-((y[i]-y[i-1])/h[i-1])); } d[m-1]=6/h[m-2]*(y3-(y[m-1]-y[m-2])/h[m-2]); } void scyt::zgf() { u[m-2]=1; for(int i=0;i<m;i++) { b[i]=2; } c[0]=1; for(i=1;i<m-1;i++) { c[i]=q[i-1]; } //........................................ l[0]=b[0]; for(i=0;i<m-1;i++) { r[i]=c[i] / l[i]; l[i+1]=b[i+1] - (u[i] * r[i]); } o[0]=d[0] / l[0]; for(i=1;i<m;i++) { o[i]=(d[i]-u[i-1]*o[i-1]) / l[i]; } M[m-1]=o[m-1]; for(i=m-2;i>=0;i--) { M[i]=o[i]-r[i] * M[i+1]; } cout<<m<<"个点的导数分别是:"<<endl; for(i=0;i<m;i++) { cout<<"M"<<i+1<<"="; cout<<M[i]<<endl; } //M的求出。。。。。。追赶法调用完毕 } void scyt::qiujie() { float S; for(;;) { float f; cout<<"请输入待求x的(输入1000)时退出:"; cin>>f; if(f==1000) break; for(int i=0;i<m;i++) { if(f>x[i]&&f<x[i+1]) { S=pow((x[i+1]-f),3)*M[i]/(6*h[i]) + pow(f-x[i],3)*M[i+1]/(6*h[i]) + (x[i+1]-f)*(y[i]-h[i]*h[i]*M[i]/6)/h[i] + (f-x[i])*(y[i+1]-h[i]*h[i]*M[i+1]/6)/h[i]; cout<<"S["<<f<<"]="<<S<<endl; } } } } scyt::~scyt() { delete []x,y,d,h,u,q,a,b,c,l,r,o,M; }
以下是三次样条算法的C语言实现代码: ```c #include <stdio.h> #include <stdlib.h> #include <math.h> #define N 100 double x[N], y[N], h[N], b[N], u[N], v[N], z[N], c[N], d[N]; int n; void spline() { int i, k; double p, qn, sig, un; for (i = 1; i < n; i++) h[i] = x[i] - x[i - 1]; for (i = 2; i < n; i++) { sig = h[i] / (h[i] + h[i - 1]); p = sig * u[i - 1] + 2.0; u[i] = (sig - 1.0) / p; b[i] = (y[i] - y[i - 1]) / h[i] - (y[i - 1] - y[i - 2]) / h[i - 1]; b[i] = (6.0 * b[i] / (h[i] + h[i - 1]) - sig * b[i - 1]) / p; } qn = un = 0.0; c[n - 1] = 0.0; for (k = n - 2; k >= 1; k--) { c[k] = u[k] * c[k + 1] + b[k]; z[k] = (c[k + 1] - c[k]) / h[k]; qn += (y[k + 1] - y[k]) / (x[k + 1] - x[k]); un += z[k] * h[k]; } z[0] = z[1]; c[0] = 3.0 * qn / h[1] - un * h[1] / 6.0; d[0] = -c[0] / h[0]; d[n - 1] = 0.0; for (k = 1; k < n; k++) { d[k] = (c[k] - c[k - 1]) / h[k]; b[k] = h[k] * (2.0 * c[k] + c[k - 1]) / 6.0 + (y[k] - y[k - 1]) / h[k]; } } double f(double t) { int i; double p; i = 0; while (x[i + 1] < t) i++; p = (t - x[i]) / h[i]; return ((1.0 - p) * y[i] + p * y[i + 1] + p * (1.0 - p) * (c[i] * (1.0 - p) + c[i + 1] * p) * h[i] * h[i] / 6.0); } int main() { int i; double t; printf("请输入数据点个数n:"); scanf("%d", &n); printf("请输入数据点的x坐标和y坐标:\n"); for (i = 0; i < n; i++) scanf("%lf%lf", &x[i], &y[i]); spline(); printf("请输入点的x坐标:"); scanf("%lf", &t); printf("在x=%g处的函数为%g\n", t, f(t)); return 0; } ```

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值