数值分析三次样条插值
“数值分析”中的三次样条插值实践教学探讨
性、 主动性、 灵活性得不到充分发挥 ， 难以提高探索 和 获取样 条 知识 的能力 。
三是 对案 例教 学 中的数 据缺乏 背 景 了解 。虽 然 在个 别教 材 中 ， 添加 了三次样 条实 际应 用案 例 ， 但仅 靠课 堂讲 授和上 机 实验是 无法 让学 生对 案例 中 的实
“ 数值 分 析 ＂ 中的 三次 样 条 插值 实践 教 学 探 讨 木
徐 圣 兵
（ 广东工业大学应用数学学 院 ， 广东 广州 ， １０６ ５０ ０ ） 摘 要： 文章 阐述 了“ 数值 分析” 中的三次样条插值教学 中加入 实践教 学环 节的必要性 ， 通过具体 实例介绍 了作 并 者在 三次样条插值教 学中进行 了实践教 学的尝试和探 讨 ， 旨在对今后 的三次样条插 值 的教 学起到一定 的
等方面要符合培养 目标的要求 ， 要和相关课程相匹 配 。实践 教学是 专 业 教学 的重 要 组 成部 分 ， 培养 是
学生 自学 、 独立 分析 问 题 和 解 决 问题 能 力 的 主要 手 段 ， 实践性 较 强课程 的一 种重 要教 学手段 。 是
环 节 已成为 当务 之急 。
一
二 、实 践 教 学
针 对 上述 问题 ， 了培 养学 生 自主分 析 问题 、 为 解
、
三次样 条插值教 学 中存在 的 问题
决 问题 的 能力 ， 发 学生学 习三次样 条插 值 的兴趣 ， 激 加 深他 们对 三次 样 条 知识 的理 解 和认 识 ， 文提 出 本
了在三 次样 条教 学 中加 入 实践教 学环 节 的想法 。 实 践教 学 的含 义包 括 实验 、 习 、 会 实 践 、 实 社 课
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;}。
- 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