数值分析三次样条插值
“数值分析”中的三次样条插值实践教学探讨
性、 主动性、 灵活性得不到充分发挥 ， 难以提高探索 和 获取样 条 知识 的能力 。
三是 对案 例教 学 中的数 据缺乏 背 景 了解 。虽 然 在个 别教 材 中 ， 添加 了三次样 条实 际应 用案 例 ， 但仅 靠课 堂讲 授和上 机 实验是 无法 让学 生对 案例 中 的实
“ 数值 分 析 ＂ 中的 三次 样 条 插值 实践 教 学 探 讨 木
徐 圣 兵
（ 广东工业大学应用数学学 院 ， 广东 广州 ， １０６ ５０ ０ ） 摘 要： 文章 阐述 了“ 数值 分析” 中的三次样条插值教学 中加入 实践教 学环 节的必要性 ， 通过具体 实例介绍 了作 并 者在 三次样条插值教 学中进行 了实践教 学的尝试和探 讨 ， 旨在对今后 的三次样条插 值 的教 学起到一定 的
等方面要符合培养 目标的要求 ， 要和相关课程相匹 配 。实践 教学是 专 业 教学 的重 要 组 成部 分 ， 培养 是
学生 自学 、 独立 分析 问 题 和 解 决 问题 能 力 的 主要 手 段 ， 实践性 较 强课程 的一 种重 要教 学手段 。 是
环 节 已成为 当务 之急 。
一
二 、实 践 教 学
针 对 上述 问题 ， 了培 养学 生 自主分 析 问题 、 为 解
、
三次样 条插值教 学 中存在 的 问题
决 问题 的 能力 ， 发 学生学 习三次样 条插 值 的兴趣 ， 激 加 深他 们对 三次 样 条 知识 的理 解 和认 识 ， 文提 出 本
了在三 次样 条教 学 中加 入 实践教 学环 节 的想法 。 实 践教 学 的含 义包 括 实验 、 习 、 会 实 践 、 实 社 课
C _数值分析_三次样条插值_自动选取步长梯形法_ROMBERG求积法_列主元高斯消去法_列主元LU分解法_JACOBI迭
//系数矩阵 //右端项 //中间项 //输出 //选取列主元的比较器
int i,j,k;
//计数器
void main() {
cout << "请输入线性方程组（ai1，ai2，ai3......ain, yi）："<<endl; for ( i = 0; i < N ;i++) {
for (int j = 0; j< N ;j++ ) cin >> A[i][j];
A[i][j] = A[i][j] - T * A[k][j]; } } } X[N-1] = B[N-1]/A[N-1][N-1]; for (i = N-2; i >=0 ; i--) {
6
double Temp = 0; for (int j = i+1; j<N ;j++)
Temp = Temp + A[i][j] * X[j]; X[i] = (B[i] - Temp) /A[i][i]; } cout << "线性方程组的解（X1，X2，X3......Xn）为："<<endl; for( i = 0; i < N ;i++) { cout << X[i] <<" "; } } 运行结果截图：
double fun(double a) {
return 2/( 1+a*a ); } double SelfSelLength(double R_a,double R_b,double e) {
double h = (R_b-R_a)/2; double R1 = (fun(R_a)+fun(R_b)) * h; int n = 1; double R0; double S; double E; do //每当误差值不符合要求时，计算下一个 result 值 {
(完整)三次样条插值的C程序(很全啊)
三次样条插值 C/C++程序（自己整理的) 具体推导看书<<数值分析〉〉 code:#include <iostream> using namespace std;(完整)三次样条插值的 C 程序(很全啊)const int MAXN = 100;int n； double x[MAXN］, y[MAXN］; //下标从 0。
n double alph［MAXN], beta[MAXN]， a[MAXN］， b[MAXN］； double h［MAXN]; double m［MAXN]; //各点的一阶导数；inline double sqr（double pa） { return pa * pa;｝double sunc(double p, int i) {［i]return (1 + 2 ＊ (p — x[i］) / (x[i + 1] - x［i］)） ＊ sqr((p - x[i + 1]) / (x[i + 1] — x［i］）） ＊ y+ （1 + 2 ＊ (p — x[i + 1]) / （x［i］ - x［i + 1]）） * sqr（(p - x［i]) / (x［i + 1] — x ［i］)) * y［i + 1]+ (p - x［i］) ＊ sqr（(p — x[i + 1］） / (x［i] - x［i + 1］)) ＊ m［i] + (p — x[i + 1]） ＊ sqr（(p — x［i]) / (x［i + 1］ — x［i])) * m[i + 1]； }int main() { int i, j；double xx;(完整)三次样条插值的 C 程序(很全啊)freopen（"threeInsert.in"， "r"， stdin)；scanf(”％d"， &n）；for （i = 0； i <= n； i++） scanf（”%lf％lf”, &x［i］, ＆y[i］）；// scanf（"％lf%lf”, ＆m［0］, &m[n］)；for （i = 0; i <= n - 1； i++) h[i］ = x[i + 1］ - x[i];//第一种边界条件//alph［0］ = 0； alph［n］ = 1； beta［0］ = 2 * m［0]； beta［n］ = 2 ＊ m[n］；//第二种边界条件alph［0］ = 1； alph［n] = 0； beta[0] = 3 * （y[1] — y［0］） / h［0]; beta［n］ = 3 * （y［n] - y[n — 1］ / h ［n — 1]）;for （i = 1； i 〈= n - 1； i++） ｛ alph[i］ = h［i — 1] / (h［i - 1] + h[i]）； beta［i] = 3 * （（1 - alph[i］) * (y［i］ — y［i — 1］） / h[i - 1］ + alph［i］ * （y[i + 1］ — y［i］） / h［i])；｝ a［0］ = — alph［0] / 2; b[0］ = beta[0］ / 2;for (i = 1； i <= n； i++) { a[i］ = — alph［i］ / （2 + （1 — alph[i］) * a[i — 1］）; b［i] = （beta［i] - （1 — alph[i］) * b［i - 1］) / (2 + （1 - alph［i]) ＊ a[i — 1］);} m［n + 1］ = 0；for (i = n; i 〉= 0; i--) ｛ m[i] = a［i］ * m[i + 1］ + b［i］;} scanf（"％lf”, &xx)；for (i = 0; i < n； i++) { if （xx 〉= x［i] ＆＆ xx <= x[i + 1］) break;｝ printf(”%lf\n", sunc(xx， i));return 0； }(完整)三次样条插值的 C 程序(很全啊)#include<iostream〉 ＃include<iomanip〉 using namespace std；const int MAX = 50; float x［MAX］, y[MAX］, h[MAX]；//变量设置：x 为各点横坐标；y 为各点纵坐标；h 为步长 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++）(完整)三次样条插值的 C 程序(很全啊)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〈〈"请输入已知断点个数：”;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;}(完整)三次样条插值的 C 程序(很全啊)for(i = 0； i < n； i++) //求步长；其数组值较之 x，y 个数少一 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"；//显示数据为 Y0' 至 Yn’，即断点的一阶导数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”;//显示数据为 Y0” 至 Yn”，即断点的二阶导数 cin>>f0>>f1; c[0] = a［n］ = 0； fxym[0] = 2*f0; fxym[n] = 2*f1; break； default:cout<<"不可用\n"；//待定｝;//switch(完整)三次样条插值的 C 程序(很全啊)for（i = 1; i 〈 n; i++)fxym[i］ = 6 * f（i—1， i， i+1）；//调用差分函数（only！）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）；//调用弯矩函数（only！）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）;//通过操作器 setprecision()设置有效位数；其为头文件 〈iomanip。
三次样条插值的方法和思路
三次样条插值的方法和思路摘要：1.三次样条插值的基本概念2.三次样条插值的数学原理3.三次样条插值的实现步骤4.三次样条插值的优缺点5.三次样条插值在实际应用中的案例正文：在日常的科学研究和工程应用中，我们经常会遇到需要对一组数据进行插值的问题。
插值方法有很多，其中三次样条插值是一种常见且有效的方法。
本文将从基本概念、数学原理、实现步骤、优缺点以及实际应用案例等方面，全面介绍三次样条插值的方法和思路。
一、三次样条插值的基本概念三次样条插值（Cubic Spline Interpolation）是一种基于分段多项式的插值方法。
它通过在各个节点上构建一条三次多项式曲线，使得这条曲线在节点之间满足插值条件，从而达到拟合数据的目的。
二、三次样条插值的数学原理三次样条插值的数学原理可以分为两个部分：一是分段三次多项式的构建，二是插值条件的满足。
1.分段三次多项式的构建假设有一组数据点序列为（x0，y0），（x1，y1），（x2，y2），（x3，y3），我们可以将这些数据点连接起来，构建一条分段三次多项式曲线。
分段三次多项式在每个子区间上都是一个三次多项式，它们之间通过节点值进行连接。
2.插值条件的满足为了使分段三次多项式在节点之间满足插值条件，我们需要在每个子区间上满足以下四个条件：（1）端点条件：三次多项式在区间的端点上分别等于节点值；（2）二阶导数条件：三次多项式在区间内的二阶导数等于节点间的斜率；（3）三阶导数条件：三次多项式在区间内的三阶导数等于节点间的曲率；（4）内部点条件：三次多项式在区间内部满足插值函数的连续性。
通过求解这四个条件，我们可以得到分段三次多项式的系数，从而实现插值。
三、三次样条插值的实现步骤1.确定插值节点：根据数据点的位置，选取合适的节点；2.构建分段三次多项式：根据节点值和插值条件，求解分段三次多项式的系数；3.计算插值结果：将待插值点的横坐标代入分段三次多项式，得到插值结果。
第4章 插值法(3)
第4章 插值逼近
相应也要求样条插值函数s(x)也具有周期性,故在端 点要求满足条件
《 数 值 分 析 》
s′( x0 ) = s′( xn ), s′′( x0 ) = s′′( xn )
(4―50)
第4章 插值逼近
hk uk = hk −1 + hk
《 数 值 分 析 》
6 yk +1 + yk yk − yk −1 λk = ( − ) hk −1 + hk hk hk −1
那么
hk −1 1 − uk = hk −1 + hk
(1 − uk )mk −1 + 2mk Байду номын сангаас uk mk +1 = λk k = 1, 2,L , n − 1
(4―44)
第4章 插值逼近
从而解出Ak和Bk,即
yk +1 − yk hk Ak = − ( mk +1 − mk ) hk 6
《 数 值 分 析 》
(4―45)
hk2 Bk = yk − mk 6
(4―46)
由式(4―43)可看出三次样条插值函数s(x)仅与 mk、m k+1有关系,因此只要求得各个mk,则各个子区 间［xk,x k+1］上的三次样条函数也就确定了。下面 介绍求mk的方法。
① 函数y=f(x)在两端点x0及xn处的导数y′0和y′n为已 知。此时要求
′ s′( x0 ) = y0 ,
′ s′( xn ) = yn
由式(4―48)和(4―49)得到
数值计算方法( 三次样条插值)
u xj hj
分段三次Hermite插值算法
则 v A1 y j 1 A2 y j B1 f j1 B2 f j
算法： 1.输入x j , f j , f j （j 0,1,...,n）； 2.计算插值 (1)输入插值点u; (2)对于j 1,2,...,n做 如果u x j 则计算A1 , A2 , B1 , B2 ; v A1 f j 1 A2 f j B1 f j1 B2 f j; 3.输出u , v。
三次样条插值
于是由Taylor展示有 s( x) s( xi ) s( xi )(x xi ) s( xi ) s( xi ) 2 ( x xi ) ( x xi )3 2! 3! M M Mi yi s( xi )(x x j ) i ( x xi ) 2 i 1 ( x xi )3 2! 3!( xi 1 xi )
2M 0 M 1 6 f [ x0 , x0 , x1 ]
三次样条插值
同理(2)式中令i n得 M n 1 2M n 6 f [ xn 1 , xn , xn ] 即有 2M 0 M 1 6 f [ x0 , x0 , x1 ] ） i M i 1 2M i i M i 1 6 f [ xi 1 , xi , xi 1 ] （i 1,2,...,n 1 M 2M 6 f [ x , x , x ] n n 1 n n n 1
三次样条插值
对于待定系数a j , b j , c j .d j j 1,2,...n,即4n个未知系数，
而插值条件为 n 2个，还缺两个，因此须 4 给出两个 条件称为边界条件，有 以下三类： 第一类 已知两端点的一阶导数 s( x0 ) f ( x0 ) m0 s( xn ) f ( xn ) mn
三次样条插值方法的应用
CENTRAL SOUTH UNIVERSITY数值分析实验报告三次样条插值方法的应用一、问题背景分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。
样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。
下面我们讨论最常用的三次样条函数及其应用。
二、数学模型样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。
设区间[]b ,a 上给定有关划分b x x n =<<<=Λ10x a ，S 为[]b ,a 上满足下面条件的函数。
● )(b a C S ,2∈；● S 在每个子区间[]1,+i i x x 上是三次多项式。
则称S 为关于划分的三次样条函数。
常用的三次样条函数的边界条件有三种类型：● Ⅰ型 ()()n n n f x S f x S ''0'',==。
● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。
● Ⅲ型 ()()Λ3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。
鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。
三、算法及流程按照传统的编程方法，可将公式直接转换为MATLAB可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB在矩阵运算上的优势。
三次样条插值计算算法
/* 三次样条插值计算算法*/#include "math.h "#include "stdio.h "#include "stdlib.h "/*N：已知节点数N+1R：欲求插值点数R+1x，y为给定函数f(x)的节点值{x(i)} (x(i) <x(i+1)) ,以及相应的函数值{f(i)} 0 <=i <=NP0=f(x0)的二阶导数；Pn=f(xn)的二阶导数u:存插值点{u(i)} 0 <=i <=R求得的结果s(ui)放入s[R+1] 0 <=i <=R返回0表示成功，1表示失败*/int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[]){/*声明局部变量*/double *h; /*存放步长:{hi} 0 <=i <=N-1 */double *a; /*存放系数矩阵{ai} 1 <=i <=N ；分量0没有利用*/ double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0 <=i <=N-1 */double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0 <=i <=N */double *af; /*存放系数矩阵{a(f)i} 1 <=i <=N ；*/double *ba; /*存放中间结果0 <=i <=N-1*/double *m; /*存放方程组的解{m(i)} 0 <=i <=N ；*/int i,k;double p1,p2,p3,p4;/*分配空间*/if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);/*第一步：计算方程组的系数*/for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++)a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++)c[k]=1-a[k];for(k=1;k <N;k++)g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]); c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;/*第二步：用追赶法解方程组求{m(i)} */ba[0]=c[0]/2;g[0]=g[0]/2;for(i=1;i <N;i++){af[i]=2-a[i]*ba[i-1];g[i]=(g[i]-a[i]*g[i-1])/af[i];ba[i]=c[i]/af[i];}af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N]; /*P110 公式：6.32*/ for(i=N-1;i> =0;i--)m[i]=g[i]-ba[i]*m[i+1];/*第三步：求值*/for(i=0;i <=R;i++){/*判断u(i)属于哪一个子区间，即确定k */if(u[i] <x[0] || u[i]> x[N]){/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 1;}k=0;while(u[i]> x[k+1])k++;//p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); //p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}/*释放空间*/free(h);free(a);free(c);free(g);free(af);free(ba);free(m);return 0;}void main(){int N,R;double *x,*y,*u,*s;double P0,Pn;int i;/*验证算法：*/N=7;R=6;/*分配空间*/if(!(x=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(y=(double*)malloc((N+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(u=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}if(!(s=(double*)malloc((R+1)*sizeof(double)))){printf( "malloc error!\n ");exit(1);}x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9 463;u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;P0=-0.4794;Pn=-0.9463;if(!SPL( N, R, x, y, P0, Pn, u, s)){/*打印结果*/printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\ns= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf( "\nsin= ");for(i=0;i <=R;i++)printf( "%9.5f ",sin(u[i]));}/*释放空间*/free(x);free(y);free(u);free(s);}/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/ //输入327.7 4.128 4.329 4.130 3.013.0 -4.0//输出输出三次样条插值函数：1: [27.7 , 28]13.07*(x - 28)^3 + 0.22*(x - 27.7)^3+ 14.84*(28 - x) + 14.31*(x - 27.7)2: [28 , 29]0.066*(29 - x)^3 + 0.1383*(x - 28)^3+ 4.234*(29 - x) + 3.962*(x - 28)3: [29 , 30]0.1383*(30 - x)^3 - 1.519*(x - 29)^3+ 3.962*(30 - x) + 4.519*(x - 29)//三次样条插值函数#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";//待定};//switchfor(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;}。
研究生数值分析 样条插值
x
-0.46 -0.40 -0.36 -0.30 -0.26 -0.20 -0.16 -0.10 -0.06 -0.00
1 1 25x2
L10 (x)
0.15898 0.24145 0.20000 0.19999 0.23585 0.18878 0.30769 0.23535
0.37175 0.31650
f (x)
L1(x)
yi 1
x xi xi1 xi
yi
x xi1 xi xi1
这种分段低次插值称为分段线性插值。
在几何上就是用折线代替曲线，故分段线 性插值又称为折线插值。
类似地，为求f(x)的近似值，也可选取距点x
最近的3个节点 xi1, xi , xi1 进行二次插值，即取
f
(x)
L2 (x)
i 1
[ yk
k i1
i 1
(
j i 1
x xj xk x j
)]
jk
这种分段低次插值叫分段二次插值。
在几何上就是用分段抛物线代替曲线，故分 段二次插值又称为抛物线插值。
3、三次样条插值 对于给定的n+1个节点，求函数的近似值，可以
作 n次插值多项式，当n较大时，高次插值不仅计算 复杂，而且还可能出现高阶导数不一致收敛的现象;
若采用分段插值，虽计算简单，也具有一致收 敛性，但光滑性比较差.
有些实际问题，比如：船体放样，飞机的机翼 设计等要求二阶光滑度(有二阶的连续导数)。过去， 工程师制图时，往往用一根富有弹性的木条（称为 样条），把它用压铁固定在样点上，其他地方让它 自由弯曲，然后画一条曲线，称为样条曲线。
它实际上是由分段三次曲线连接而成，在连接 点处有二阶连续导数。我们对工程师描绘的样条曲 线，抽象成数学模型，得出的函数称为样条函数， 它实质上是分段多项式的光滑连接。
数值分析课程设计--三次样条插值
《数值分析》课程设计三次样条插值算法院（系）名称信息工程学院专业班级 09普本信计1班学号 090111073学生姓名宣章然指导教师孔繁民2012年06月08日数值分析课程设计评阅书课程设计任务书2008—2009学年第二学期专业班级： 09普本信计1班学号： 060111060 姓名：宣章然课程设计名称：数值分析设计题目：三次样条插值完成期限：自 2012 年 6 月 8 日至 2012 年 6 月 13 日共 1 周设计依据、要求及主要内容：一、设计目的熟练掌握三次样条插值算法的原理和推导过程，并且能够应用Matlab软件编写相应的程序和使用Matlab软件函数库软件。
二、设计要求(1)用Matlab函数库中相应函数对选定的问题，求出具有一定精度的结果。
(2)使用所用的方法编写Matlab程序求解,对数值结果进行分析。
(3)对于使用多个方法解同一问题的，在界面上设计成菜单形式。
三、设计内容首先构造三次样条插值函数的定义和一般特征,并对实例问题进行实例分析,并总结四、参考文献[1] 黄明游,冯果忱.数值分析[M].北京:高等教育出版社,2008.[2] 马东升,雷勇军.数值计算方法[M].北京:机械工业出版社,2006.[3] 石博强,赵金.MATLAB数学计算与工程分析范例教程[M].北京:中国铁道出版社.2005.[4]郝红伟，MATLAB 6，北京，中国电力出版社，2001[5]姜健飞，胡良剑，数值分析及其MATLAB实验，科学出版社，2004[6]薛毅，数值分析实验，北京工业大学出版社，2005 计划答辩时间：2012年6月18日指导教师（签字）：教研室主任（签字）：批准日期：年月三次样条插值摘 要分段低次样条插值虽然计算简单、稳定性好、收敛性有保证且易在电子计算机上实现，但只能保证各小段曲线在连接处的连续性，不能保证整件曲线的光滑性。
利用样条插值，既可保持分段低次插值多项式，又可提高插值函数光滑性。
三次样条插值算法详解
三次样条插值算法要求数据点数量较多，且在某些情况下可能存在数值不稳定性，如数据 点过多或数据点分布不均等情况。此外，该算法对于离散数据点的拟合效果可能不如其他 插值方法。
对未来研究的展望
01
02
03
改进算法稳定性
针对数值不稳定性问题， 未来研究可以探索改进算 法的数值稳定性，提高算 法的鲁棒性。
3
数据转换
对数据进行必要的转换，如标准化、归一化等， 以适应算法需求。
构建插值函数
确定插值节点
根据数据点确定插值节点，确保插值函数在节点处连续且光滑。
构造插值多项式
根据节点和数据点，构造三次多项式作为插值函数。
确定边界条件
根据实际情况确定插值函数的边界条件，如周期性、对称性等。
求解插值函数
求解线性方程组
06
结论
三次样条插值算法总结
适用性
三次样条插值算法适用于各种连续、光滑、可微的分段函数插值问题，尤其在处理具有复 杂变化趋势的数据时表现出色。
优点
该算法能够保证插值函数在分段连接处连续且具有二阶导数，从而在插值过程中保持数据 的平滑性和连续性。此外，三次样条插值算法具有简单、易实现的特点，且计算效率较高 。
根据数据点的数量和分布，合理分段，确保 拟合的精度和连续性。
求解线性方程组
使用高效的方法求解线性方程组，如高斯消 元法或迭代法。
结果输出
输出拟合得到的插值函数，以及相关的误差 分析和图表。
03
三次样条插值算法步骤
数据准备
1 2
数据收集
收集需要插值的原始数据点，确保数据准确可靠。
数据清洗
对数据进行预处理，如去除异常值、缺失值处理 等。
数值分析实验报告-插值、三次样条
实验报告：牛顿差值多项式&三次样条问题：在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数21()25f x x作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。
实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。
应用所编程序解决实际算例。
实验要求：1． 认真分析问题，深刻理解相关理论知识并能熟练应用； 2． 编写相关程序并进行实验； 3． 调试程序，得到最终结果； 4． 分析解释实验结果； 5． 按照要求完成实验报告。
实验原理：详见《数值分析 第5版》第二章相关内容。
实验内容：（1）牛顿插值多项式1.1 当n=10时：在Matlab 下编写代码完成计算和画图。
结果如下：代码：clear allclcx1=-1:0.2:1;y1=1./(1+25.*x1.^2);n=length(x1);f=y1(:);for j=2:nfor i=n:-1:jf(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1));endendsyms F x p;F(1)=1;p(1)=y1(1);for i=2:nF(i)=F(i-1)*(x-x1(i-1));p(i)=f(i)*F(i);endsyms PP=sum(p);P10=vpa(expand(P),5);x0=-1:0.001:1;y0=subs(P,x,x0);y2=subs(1/(1+25*x^2),x,x0);plot(x0,y0,x0,y2)grid onxlabel('x')ylabel('y')P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x^10+494.91*x^8-9.5065e-14*x^7-381.43*x^6-8.504e-14*x^5+123.36* x^4+2.0202e-14*x^3-16.855*x^2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。
三次样条插值求导法
三次样条插值求导法三次样条插值法是一种常用的数值分析方法，用于近似插值实现平滑曲线的拟合。
它的优点在于可以保持原始数据的特性，同时能够降低数据间的噪声干扰，使得插值的结果更加准确。
本文将介绍三次样条插值法的原理、算法以及应用方面的指导意义。
首先，我们需要了解三次样条插值法的基本原理。
三次样条插值法通过在相邻数据点之间构造三次多项式来近似拟合原始数据。
这些三次多项式满足一定的光滑性条件，使得插值结果的曲线平滑而连续。
在三次样条插值中，每个数据点都对应一个三次多项式，并且相邻多项式之间的导数和二阶导数必须相等，以保证曲线的平滑性。
接下来，我们将介绍三次样条插值法的算法步骤。
首先，我们需要确定每个数据点对应的三次多项式。
为了满足光滑性条件，我们需要计算每个数据点处的导数值。
这可以通过求解一个线性方程组来实现，其中方程的个数等于数据点的个数。
解得导数值之后，我们就可以得到每个数据点对应的三次多项式的系数。
然后，我们需要利用这些系数来计算在数据点之间的插值结果。
为了实现这一点，我们可以利用三次多项式的性质，通过给定的数据点和对应的三次多项式系数，来计算在两个相邻数据点之间的插值结果。
最后，我们需要通过合理的选择数据点以及插值节点的间距，来获得更加准确的三次样条插值结果。
一般来说，数据点的选择应尽量满足曲线的变化趋势，以反映原始数据的特点。
此外，插值节点的间距也需要经过合理的选择，以保证插值结果的准确性。
三次样条插值法在实际应用中有着广泛的意义和指导价值。
首先，它可以用于光滑曲线的拟合，将离散的数据点进行连续化处理，使得数据的绘图和分析更加方便。
其次，它可以用于数据的插值预测，通过已有的数据点来预测未知数据点的取值。
此外，三次样条插值法还可以在数字图像处理中用于图像的平滑和插值填充，从而改善图像的质量和美观度。
综上所述，三次样条插值法是一种有效的数值分析方法，可以用于实现平滑曲线的拟合和数据的插值预测。
通过了解其原理、算法以及应用方面的指导意义，我们可以更好地理解和应用这一方法，从而提高数据处理和分析的准确性和效率。
第二章三次样条插值
mk 1 2mk
hk 1 hk hk 1
k 1
3( hk yk1 yk hk1 yk yk1 )
hk hk 1
hk
hk hk 1
hk 1
k mk1 2mk k mk 1 gk
k
hk
hk hk 1
k
hk 1 hk hk 1
gk
3(k
yk yk 1 hk 1
根据S ( x)在[a, b]上二阶导数连续 在节点xj j (1, 2, , n 1)处就应满足连续性条件
S(x j 0) S(x j 0) S ' (x j 0) S '(x j 0) S"(x j 0) S"(x j 0)
共(n+1)+(3n-3)=4n-2个条件，因此还需要 两个条件才能确定S(x)
注：三次样条与分段 Hermite 插值的根本区别在于S(x)自 身光滑，不需要知道 f 的导数值（除了在2个端点可能需 要）；而Hermite插值依赖于f 在所有插值点的导数值。
f(x)
H(x)
S(x)
要求出S(x)，则在每个小区间上 [x j , x要j1确] 定4个 待定系数，共有n个小区间，所以应确定4n个参 数。
可在区间端点a,b上各加一个条件（边界条件）， 具体要根据实际问题要求给定；
1. 已知两端的一阶导数值
S(x0 ) f0
S(xn ) fn
2. 两端的二阶导数已知
S(x0 ) f0 S"(xn ) fn
其特殊情况为
S(x0 ) S(xn ) 0
3. 当f(x)是xn x0为周期的周期函数时，则要求 S(x)也是周期函数，这时边界条件应满足：
三次样条插值c++代码实现及注释
一、引言在计算机编程和数据处理领域，插值是一种常见的数值分析方法，用于在已知数据点之间估算未知点的数值。
而三次样条插值是插值方法中的一种重要技术，它可以在使用较少插值节点的情况下，实现更为平滑和精确的插值结果。
本文将着重探讨三次样条插值的原理和C++代码实现，并给出详细的注释和解释。
二、三次样条插值的原理三次样条插值是一种分段插值方法，它将整个插值区间分割为若干个小区间，每个小区间内采用三次多项式进行插值。
这样做的好处是可以在每个小区间内实现更为细致和精确的插值，从而提高插值的准确性和平滑性。
而三次样条插值的核心在于确定每个小区间内的三次多项式的系数，一般采用自然边界条件进行求解。
在具体实现中，我们需要先对给定的插值节点进行排序，并求解出每个小区间内的三次多项式系数。
最终将这些系数整合起来，就可以得到整个插值区间的三次样条插值函数。
三、C++代码实现及注释接下来，我们将给出使用C++语言实现三次样条插值的代码，并对每个关键步骤进行详细注释和解释。
```cpp// include necessary libraries#include <iostream>#include <vector>using namespace std;// define the function for cubic spline interpolationvector<double> cubicSplineInterpolation(vector<double> x, vector<double> y) {// initialize necessary variables and containersint n = x.size();vector<double> h(n-1), alpha(n), l(n), mu(n), z(n), c(n), b(n), d(n);vector<double> interpolatedValues;// step 1: calculate the differences between x valuesfor (int i = 0; i < n-1; i++) {h[i] = x[i+1] - x[i];}// step 2: calculate alpha valuesfor (int i = 1; i < n-1; i++) {alpha[i] = (3/h[i]) * (y[i+1] - y[i]) - (3/h[i-1]) * (y[i] - y[i-1]); }// step 3: calculate l, mu, and z valuesl[0] = 1;mu[0] = 0;z[0] = 0;for (int i = 1; i < n-1; i++) {l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*mu[i-1];mu[i] = h[i]/l[i];z[i] = (alpha[i] - h[i-1]*z[i-1])/l[i];}l[n-1] = 1;z[n-1] = 0;c[n-1] = 0;// step 4: calculate coefficients for the cubic polynomials for (int j = n-2; j >= 0; j--) {c[j] = z[j] - mu[j]*c[j+1];b[j] = (y[j+1] - y[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3;d[j] = (c[j+1] - c[j])/(3*h[j]);}// step 5: interpolate values using the cubic polynomials for (int i = 0; i < n-1; i++) {double xi = x[i];while (xi < x[i+1]) {double dx = xi - x[i];double interpolatedValue = y[i] + b[i]*dx + c[i]*dx*dx + d[i]*dx*dx*dx;interpolatedValues.push_back(interpolatedValue);xi += 0.1; // adjust the step size for finer interpolation }}return interpolatedValues;}// main function for testing the cubic spline interpolation int main() {vector<double> x = {1, 2, 3, 4, 5};vector<double> y = {3, 6, 8, 10, 15};vector<double> interpolatedValues = cubicSplineInterpolation(x, y);for (int i = 0; i < interpolatedValues.size(); i++) {cout << "Interpolated value " << i << " : " << interpolatedValues[i] << endl;}return 0;}```四、总结与展望通过本文的学习，我们了解了三次样条插值的原理和C++代码实现。
计算方法大作业——三次样条插值
计算方法上机报告
此完成所有数据的输入。继续按 Enter 键会出现提示“选择封闭方程组的边界条件: 第 一类边界条件输入 1，第二类边界条件输入 2，第三类边界条件输入 3。 ”根据已知情况 选择相应的边界条件，若为自然三次样条插值，则选 1，并将插值区间两端点的二阶导 数值设置为 0。输入完成之后按 Enter 开始求解，程序运行结束后命令窗口会显示要求 的三次样条插值函数，同时会出现该插值函数以及插值节点的图像，便于直接观察。 2.3 算例及计算结果 (1) 《数值分析》课本第 137 页的例题 4.6.1，已知函数 y=f(x)的数值如下表，求它 的自然三次样条插值函数。 xi yi -3 7 -1 11 0 26 3 56 4 29
(2) 给定函数 f ( x)
3 x 1 1 x 0 0 x3 3 x 4
1 (1 x 1) 。取等距节点，构造牛顿插值多项式 N5(x) 1 25x 2 和 N10(x)及三次样条插值函数 S10(x)。分别将三种插值多项式与 f(x)的曲线画在同一个
N10 x
22757 10 5444 8 20216 6 17147 4 3725 2 x x x x x 1 103 11 53 139 221
将牛顿插值多项式 N5(x)和 N10(x)及三次样条插值函数 S10(x)分别与 f(x)的曲线画在 同一个坐标系上进行比较，如图 12。可以看出三次样条函数与原函数符合的非常好， 对于低次的牛顿插值多项式，与原函数的大致趋势相同，而高次的牛顿插值多项式由 于龙格现象的出现，与原函数之间相差比较大。
S ( xi ) S ( xi ), ( xi ) S ( xi ), S S ( x ) S ( x ), i i i 1, 2, , n 1
数值分析(15)样条插值
数值分析
同理，在[ xi1, xi ]也可以得到
S
" i
1
(
x
)
6x
2 xi1 h2
i 1
4xi
mi 1
6x
4 xi1 h2
i 1
2xi
mi
6( xi1 xi 2
在内节点x（i ih3i
1
1,
2,
x
)
,
( yi
n-
yi1 )
1）上，由S
" i
(
xi
)
S
" i
1
(
xi
)
1 hi 1
mi 1
化简后得到三弯矩方程
hi1Mi1 （2 hi1+hi）Mi hi Mi1 6( f xi , xi1 f xi1 , xi ) gi
(i 1, 2, , n 1)
h0 2(h0 h1 )
h1
h1
2(h1 h2 ) h2
M0 g1
M1
g2
hn2
2(hn2 hn1 )
是三次多项式; 则称S( x)为三次样条函数。x1, ..., xn1称为内节点, x0 , xn称为外节点.
数值分析
数值分析
插值条件分析
由（3）S( x)在每个[ xi , xi1]上表达式不同，故应分段构造：
S0( x)
S
(
x)
S1( x)
Sn1( x)
x [ x0 , x1]; x [ x1, x2 ];
数值分析
第五节 样条插值
样条是绘图员用于描绘光滑曲线的一种机 械器件,它是一些易弯曲材料制成的窄条或棒条. 在绘制需要通过某点的光滑曲线时,对它在这些 点的位置上“压铁”,它就被强制通过或接近图 表上确定的描绘点.“样条函数”这个术语意在 点出这种函数的图象与机械样条画出的曲线很 象.
数值分析学习课件
三次样条插值
例：已知函数y=f(x)的数表如下表所示。 已知函数y=f(x)的数表如下表所示。 y=f(x)的数表如下表所示 x f(x)
0 1
0.15
0.30
0.45
0.60
0.97800 0.91743 0.83160 0.73529
求满足边界条件
s′(0) = 0, s′(0.60) = −0.64879
已知端点二阶导数
s′′( x0 ) = f ′′( x0 ) = M 0 s′′( xn ) = f ′′( xn ) = M n
当M 0 = M n = 0 为自然边界条件 已知周期边界条件
s ( x0 ) = s ( xn ) s′( x0 + 0) = s′( x0 − 0) s′′( x = 0) = s′′( x − 0) n n
则称 s ( x ) 为区间[a, b] 对应于划分∆ 的三次样条函数。 的三次样条函数。
三次样条插值
设三次样条函数s(x) 在每个子区间[xj−1, xj ]上有表达式
s(x) = sj (x) = aj x3 + bj x2 + cj x + d j x ∈(xj−1, xj ), j =1,2...n
三次样条插值
利用三弯矩阵构造三次样条插值函数 令 S ′′( xi ) = M i (i = 0,1, 2,...n) 因为S ( x) 在[ xi , xi +1 ] 在 上是三次多项式， 上是三次多项式，所以 S ′′( x)
xi +1 − x x − xi + M i +1 hi hi
[ xi , xi +1 ] 上是线性函数， 上是线性函数，故有
数值分析作业-三次样条插值
数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。
实验函数:dt ex f xt ⎰∞--=2221)(π实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。
实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：kx012345678910 ky0.00.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29ky'0.80.2算法描述：拉格朗日插值：其中是拉格朗日基函数，其表达式为：()∏≠=--=nijj jiji xxxxxl)()(牛顿插值：))...()(](,...,,[....))(](,,[)0](,[)()(11211211----++--+-+=nnnxxxxxxxxxxfxxxxxxxfxxxxfxfxN其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[11211xxxxxfxxxfxxxfxxxxfxxfxxxfxxxfxfxxfnnnnikjikjkjijijiji三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[xi-1,xi]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131iiiiiiiiiiiiiiiiiiiiixxxhyMhMhhyxMMhhyyhxxMihxxMxS-------∈-+-+---+-+-=式中Mi=)(ixS''.因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：令⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=---+=+=+=+--++++++],,[6)(6111111111i i i i i i i i i i i i i i i i i i i ix x x f h y y h y y h h d h h h h h h λμ则Mi 满足如下n-1个方程：1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(000 （2） 给定区间两端点的二阶导数M0，Mn,即n n n M y x S M y x S =''=''=''='')(,)(000 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S （x ）也为周期函数，对S （x ）加上周期条件2,1,0),0()0()(0)(=-=+p x S x S n p p对于第一类边界条件有⎪⎪⎩⎪⎪⎨⎧--=+--=+--)(62)(6211001110n n n n n n i h y y mn h M M m h y y h M M对于第二类边界条件有⎩⎨⎧=+=+-n n n n d M M d M M 221100μλ其中n n n n nnn M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(6100001010-+-=-+-=-μλλ那么解就可以为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得0010012d M M M n =-++μλ，其中]),1[],[(6,,101010110n n nn n n x x f x x f h h d h h h h h h --+=+=+=μλ，那么解就可以为： ⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,..22n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标 %xi 插值点处的横坐标%f 求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1l=l.*(xi-X(j))/(X(i)-X(j)); end ; for j=i+1:nl=l.*(xi-X(j))/(X(i)-X(j)); end ;%拉格朗日基函数 f=f+l*Y(i); endfprintf('%d\n',f) return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi) %X 为已知数据的横坐标 %Y 为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值n=length(X);newt=[X',Y'];%计算差商表for j=2:nfor i=n:-1:1if i>=jY(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1));else Y(i)=0;endendnewt=[newt,Y'];end%计算牛顿插值f=newt(1,2);for i=2:nz=1;for k=1:i-1z=(xi-X(k))*z;endf=f+newt(i-1,i)*z;endfprintf('%d\n',f)return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi)% X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))...+M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));%三样条插值函数表达式endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6 *h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return4 三次样条插值第二类边界条件Threch2.mfunction [Sx]=Threch2(X,Y,d2y0,d2yn,xi)X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值%d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=2*d2y0;d(n+1)=2*d2yn;%赋初值A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=0;A(n+1,n)=0;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;syms x;for i=1:nSx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3);digits(4);Sx(i)=vpa(Sx(i));endfor i=1:ndisp('S(x)=');fprintf('%s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endfor i=1:nif xi>=X(i)&&xi<=X(i+1)S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i)) /(6*h(i))*(xi-X(i))^3;endenddisp('xi S');fprintf('%d,%d\n',xi,S);return5插值节点处的插值结果main3.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果');newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=[0.0,0.1,0.2,0.3,0.4];Y=[0.5000,0.5398,0.5793,0.6179,0.7554];a=linspace(0,0.4,21);NUM=21;L=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUMxi=a(i);L(i)=lang(X,Y,xi);% 拉格朗日插值N(i)=newton(X,Y,xi);%牛顿插值B(i)=normcdf(xi,0,1);%原函数S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'--r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.^2);for i=1:3 %节点数量变化次数N=2*N;t=linspace(-1,1,N+1);%插值节点ft=1./(1+25*t.^2);%插值节点函数值val=linspace(-1,1,101);for j=1:101L(j)=lang(t,ft,val(j));S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str);legend('原函数','拉格朗日插值','三次样条插值');%显示图例hold offfigureend8车门曲线main5.mclearclcX=[0,1,2,3,4,5,6,7,8,9,10];Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29]; dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i);f1(i)=(Y(i+1)-Y(i))/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);endd(1)=6*(f1(1)-dy0)/h(1);d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1);B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1B(i)=h(i)/(h(i)+h(i+1));C(i)=1-B(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1A(i,i)=2;endfor i=2:nA(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=A\d;x=zeros(1,n);S=zeros(1,n);for i=1:nx(i)=X(i)+0.5;S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)-M (i))/(6*h(i))*(x(i)-X(i))^3;endplot(X,Y,'k'); hold on;plot(x,S,'o');title('三次样条插值效果图');legend('已知插值节点','三次样条插值');hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一张图上4.5.1增加插值节点观察误差变化从上面三张图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线（注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
0
2
1
n1
1
n2
2 n1
M d 0
MM dd n2 M d 2
1 1 2 2
n1 n1 n n
di f xi2, xi1, xi
华长生制作
7
2、 三弯矩构造法
三次样条插值函数 S( x) 可以有多种表达式，有时用二阶导数
值S( xi) Mi (i 0,1,, n)
Mi
xi
表示时，使用更方便。 在力学上解释
为细M梁i 在 S处( x的) 弯矩，并且得到的弯矩与相邻两个弯矩有关，故
称用由于表S(示x)在区间的算[x法i , x为i三1](弯i 矩0,算1,法,。n 1) 上是三次多项式，
hn
n1 3
Mn
f
x0 , x1 f
xn1, xn
其中
0
h1 h1h n
1
0 ,
hn , 0 hnh0
d1
6(
f
[
x
,
0
x1]
f
x[ , n1
x
n])(h1
h
n)
1
.
可解出 M i (i 0,1,, n) ，方程组的矩阵形式为
2
hi
min hi
,M4
max x[a,b]
f (4) (x)
1in
华长生制作
16
精品课件！
精品课件！
可见S(x), S(x)和S(x)在[a,b]上一致收敛到f (x), f (x)和f (x)
显然若f (x) C4[a,b],则有 max| f (k)(x) S(k)(x)| o(h4k )
华长生制作
2
(1) S(x), S(x), S(x)都在区间[a,b]上连续 ,即S(x) C2[a,b] (2) S(x)在每个小区间[xk , xk 1 ]上都是三次多项式
则称 S(x)为区间[a,b]上的三次样条函数
(3) 如果函数 f (x)在节点x0 , x1 ,, xn处的函数值为 f (x j ) y j , j 0,1,, n
个方程.可由前面给出的任一种边界条件补充两个方程。
对于转角边界条件，可以导出两个方程
2M 0
M
n1
M1 2M n
6 (f
h1
6
hn
[ (
x f
x, ]
01
f[ n
f x
) 0 x, n1
]).
n
华长生制作
10
这样可解出 M i (i 0,1,, n) ，从而得 S( x) 的表达式,若令
5
Sk (x)是[xk , xk 1 ]上的三次样条插值多项 式,应有4个待定的系数 即要确定 S(x)必须确定 4n个待定的系数 少两个条件 并且我们不能只对插值函数在中间节点的状态进行限制
也要对插值多项式在两端点的状态加以要求
也就是所谓的边界条件:
第一类(转角)边界条件: S(x0 0) m0 S(xn 0) mn------(5)
d 0
6
h1
(
f
[
x
,
0
x
]
1
f
), 0
n
1,
0
则方程组可以写成矩阵形式
dn
6 (
hn
f
n
x x f [ , ]),
n1
n
2
0
1
2
1
n1
2
n
MMM ddd n1 M d 2
背景
L-插值（牛顿插值） 高次插值出现龙格现象
Hermite插值
分段插值 但分段线性插值在节点处不一定光滑
分段Hermite插值 但导数值不容易提取（找到）
三次样条插值（先由函数值确定导数值，再由分段 Hermite插值解决问题） 举例：
1 汽车、船的外形设计，流体力学等要求流线型（光滑）；
2 木样条的来源。
对 S( x)求导得
S( x)
Mi
( xi 1 2hi
x)2
M
i
1
(
x x 2hi
i
)2
f [xi ,
x ] i 1
hi
6
M (
i 1
M i),
由此可得
h S( 0) f [ , ] i (
x x x 2M M i
i i1 6
),
i
i1
h S( 0) f [ , ] i (
f (x j ) y j , j 0,1,, n 如果S(x)是f (x)的三次样条插值函数 ,则
S (x j ) y j , j 0, n
lim xxj
S(x)
lim
xxj
S(xj
)
lim
xx_j
S ( x)
lim
x x j
S ( x j
lim
则对任意
x 有[a,b]
f
x
sx
5 384
M 4h4
,
f
x
s x
1 8
M 4h2
,
f
x
s x
1 24
M 4h3
f
x s x
1
2
M 4h
其中hi
xi
xi1, h
max
1in
hi ,
max
1in
Sk
(x)
lim
x xk
Sk 1 ( x)
k 1,2,,n 1 k 1,2,,n 1 ------(8)
lim x xk
Sk( x)
lim
x xk
Sk1 ( x)
k 1,2,,n 1
S(x0 0) m0 S(xn 0) mn 或 S(x0 0) M 0S(xn 0) M n
华长生制作x
x
_ j
S ( x)
lim
x x j
S(x
j
)
y j , j 1, , n ) mj (未知）, j M j (未知）,
1 1,
j
,n 1,
------(2) 1
,n 1
4
共4n 2个条件
S ( x)在[a, b]上必 然是分段函数, 亦即
).
x x x M 2M i1
i i1 6
i
i1
x x 当
x [
,
i 1
]
i
时,
S( x) 的表达式可得,因此有
h S( 0) f [ , ] i1 (
).
x x x M 2M i
i1 i
6
i 1
i
利用条件 S( xi 0) S( xi 0) 得
华长生制作
9
其中
M M M d i
2
i1
i i
,i 0,1, 2, , n 1,
i1
i
h h
i1 ,
i
i
i 1 i ,
h h h h i1
i
i1
i
d x x x 6 f [ ,
].
i
i1 i, i1
上述方程组是关于 Mi 的方程组,有 n 1 个未知数,但只有 n 1
第二类(弯矩)边界条件 S(x0 0) M0 S(xn 0) M-n-----(6)
S xn 0 Sx0 0 0时，称为自然边界条件
第三类(周期)边界条件 S ( p) (x0 0) S ( p) (xn 0) ------(7)
华长生制作
S0(x)
S(x)
S1 ( x)
Sn1( x)
x [x0 , x1 ] x [x1 , x2 ]
------(3)
x [ xn1 , xn ]
Sk (x)是[xk , xk 1 ]上的(两点)三次样条插值多项式 ,满足
Sk (xj ) yj
lim x xk
华长生制作
1
什么是样条:
样条本质上是一段一段的三次多项式拼合而成的曲线 在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的
1946年,Schoenberg将样条引入数学,即所谓的样条函数
一、三次样条插值函数
定义1. a x0 , x1 ,, xn b为区间[a,b]的一个分割 如果函数 S(x)在区间[a,b]上满足条件 :
hn 1 2
fn
2
2 /3
1 2 1/3
1/3 2 1
m0 m1
2/3 m2
2 m3
g0 g1 g2 g3
解方程组得：m0
17 8
, m1
7 4
, m2
5 4
, m3