数值分析作业-三次样条插值..
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：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
三次样条插值作业题
例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表：且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
以下为Matlab 代码：%============================= % 本段代码解决作业题的例1 %============================= clear all clc% 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];LeftBoun = 0.2;RightBoun = -1;% 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1);for i = 1 : length(IndVar) - 1h(i) = IndVar(i + 1) - IndVar(i);end% 为向量μ赋值mu = zeros(1, length(h));for i = 1 : length(mu) - 1mu(i) = h(i) / (h(i) + h(i + 1));endmu(i + 1) = 1;% 为向量λ赋值lambda = zeros(1, length(h));lambda(1) = 1;for i = 2 : length(lambda)lambda(i) = h(i) / (h(i - 1) + h(i)); end% 为向量d赋值d = zeros(1, length(h) + 1);d(1) = 6 * ( (DepVar(2) - DepVar(1) ) / ( IndVar(2) - IndVar(1) ) - LeftBoun) / h(1); for i = 2 : length(h)a = ( DepVar(i) - DepVar(i - 1) ) / ( IndVar(i) - IndVar(i - 1) );b = ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) );c = (b - a) / ( IndVar(i + 1) - IndVar(i - 1) );d(i) = 6 * c;endd(i + 1) = 6 *( RightBoun - ( DepVar(i + 1) - DepVar(i) ) / ( IndVar(i + 1) - IndVar(i) ) ) / h(i);% 为矩阵A赋值% 将主对角线上的元素全部置为2A = zeros( length(d), length(d) );for i = 1 : length(d)A(i, i) = 2;end% 将向量λ的各元素赋给主对角线右侧第一条对角线for i = 1 : length(d) - 1A(i, i + 1) = lambda(i);end% 将向量d的各元素赋给主对角线左侧第一条对角线for i = 1 : length(d) - 1A(i + 1, i) = mu(i);end% 求解向量MM =A \ d';% 求解每一段曲线的函数表达式for i = 1 : length(h)Coefs_1 = M(i) / (6 * h(i));Part_1 = conv( Coefs_1, ...conv( [-1, IndVar(i + 1)], ...conv( [-1, IndVar(i + 1)], [-1, IndVar(i + 1)] ) ) ); S_1 = polyval (Part_1, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_2 = M(i + 1)/(6 * h(i));Part_2 = conv( Coefs_2, ...conv( [1, -IndVar(i)], ...conv( [1, -IndVar(i)], [1, -IndVar(i)] ) ) );S_2 = polyval (Part_2, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_3 = (DepVar(i) - M(i) * h(i)^2 / 6) / h(i);Part_3 = conv(Coefs_3, [-1, IndVar(i + 1)]);S_3 = polyval (Part_3, [IndVar(i) : 0.01 : IndVar(i + 1)]);Coefs_4 = (DepVar(i + 1) - M(i + 1) * h(i)^2 / 6) / h(i);Part_4 = conv(Coefs_4, [1, -IndVar(i)]);S_4 = polyval (Part_4, [IndVar(i) : 0.01 : IndVar(i + 1)]);S = S_1 + S_2 + S_3 + S_4;plot ([IndVar(i) : 0.01 : IndVar(i + 1)], S, 'LineWidth', 1.25)% 在样条插值曲线的相应位置标注该段曲线的函数表达式text(i - 1, polyval(Part_1, 3), ...['\itS', num2str(i), '(x)=', num2str(Coefs_1), '(', num2str( IndVar(i + 1) ), '-x)^{3}+', ...num2str(Coefs_2), '(x-', num2str( IndVar(i) ), ')^{3}+', num2str(Coefs_3), ...'(', num2str( IndVar(i + 1) ), '-x)+', num2str(Coefs_4), '(x-',num2str( IndVar(i) ), ')'], ...'FontName', 'Times New Roman', 'FontSize', 14)hold onend% 过x=1和x=2两个横轴点作垂线 %line([1, 1], [2.5, -0.5], 'LineStyle', '--');line([2, 2], [2.5, -0.5], 'LineStyle', '--');% 为x轴和y轴添加标注xlabel( '\itx', 'FontName', 'Times New Roman', ...'FontSize', 14, 'FontWeight', 'bold');ylabel( '\its(x)', 'FontName', 'Times New Roman', ...'Rotation', 0, 'FontSize', 14, 'FontWeight', 'bold');最终，三次样条插值函数s(x)表达式为：[][][]⎪⎩⎪⎨⎧∈-+-+-+--∈-+-+---∈+-++--=.3,2,)2(44.1)3(62.2)2(06.0)3(62.0,2,1,)1(62.2)2(08.0)1(62.0)2(42.0,1,0,08.0)1(06.042.0)1(06.0)(333333x x x x x x x x x x x x x x x x s曲线的图像如图所示：例2 已知函数值表：试求在区间[1,5]上满足上述函数表所给出的插值条件的三次自然样条插值函数)(x s本算法求解出的三次样条插值函数将写成三弯矩方程的形式：)()6()()6()(6)(6)(211123131j j jj j j jj j j j jj j jj x x h h M y x x h h M y x x h M x x h M x s --+--+-+-=+++++其中，方程中的系数jj h M 6，jj h M 61+，jj j j h h M y )6(2-，jjj j h h M y )6(211++-将由Matlab代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。
三次样条插值方法的应用
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在矩阵运算上的优势。
数值分析上机实验报告(插值)
数值分析第一次上机练习实验报告——Lagrange 插值与三次样条插值一、 问题的描述设()2119f x x =+, []1,1x ∈-,取15iix =-+,0,1,2,...,10i =.试求出10次Lagrange 插值多项式()10L x 和三次样条插值函数()S x (采用自然边界条件),并用图画出()f x ,()10L x ,()S x .二、 方法描述——Lagrange 插值与三次样条插值我们取15i ix =-+,0,1,2,...,10i =，通过在i x 点的函数值()2119i i f x x =+来对原函数进行插值，我们记插值函数为()g x ，要求它满足如下条件：()()21,0,1,2,...,1019i i i g x f x i x ===+ (1)我们在此处要分别通过Lagrange 插值(即多项式插值)与三次样条插值的方法对原函数()2119f x x=+进行插值，看两种方法的插值结果，并进行结果的比较。
10次的Lagrange 插值多项式为：()()10100i i i L x y l x ==∑ (2)其中：()21,0,1,2,...,1019i i iy f x i x ===+ 以及()()()()()()()()()011011......,0,1,2,...,10......i i n i i i i i i i n x x x x x x x x l x i x x x x x x x x -+-+----==----我们根据(2)进行程序的编写，我们可以通过几个循环很容易实现函数的Lagrange 插值。
理论上我们根据区间[]1,1-上给出的节点做出的插值多项式()n L x 近似于()f x ，而多项式()n L x 的次数n 越高逼近()f x 的精度就越好。
但实际上并非如此，而是对任意的插值节点，当n →+∞的时候()n L x 不一定收敛到()f x ；而是有时会在插值区间的两端点附近会出现严重的()n L x 偏离()f x 的现象，即所谓的Runge 现象。
计算方法大作业——三次样条插值
计算方法上机报告
此完成所有数据的输入。继续按 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 三次样条插值
2 三次样条插值
2.1 算法原理及程序框图 设在区间[a, b]上给定 n+1 个节点 xi（a ≤ x0 < x1 < … < xn ≤ b）,在节点 xi 处的函数 值为 yi = f(xi) (i = 0,1,…,n)。若函数 S(x)满足以下三个条件： (1) 在每个子区间[xi-1, xi] (i = 0,1,…,n)上，S(x)是三次多项式； (2) S(xi) = yi (i = 0,1,…,n)； (3) 在区间[a, b]上，S(x)的二阶导数 S”(x)连续， 则称 S(x)为函数 yi = f(x) 在区间[a, b]上的三次样条插值函数。 由定义可知 S(x)共有 4n 个待定参数，根据条件(3)可得如下 3n-3 个方程，
S x
x x i
6hi
3
M i 1
x xi 1
6hi
3
x x hi2 M i yi 1 M i 1 i 6 hi
数值计算方法三次样条插值
4.4 三次样条插值
令
A1
j1 (u )
(1
2
u
x hj
j 1
)(
u
xj hj
)2
A2
j (u )
(1
2
u
x hj
j
)(
u
x hj
j 1
)2
B1
j1 (u )
(u
u x j 1 )(
xj hj
)2
B2
j (u )
(u
x
j )(
u
x hj
j
)2
分段三次Hermite插值算法
I2(x)
I
n
(
x
)
x ( x0 , x1)
x ( x1, x2 ) ...... x ( xn1, xn )
其I中 j xxj1 xxj j yj1xxj xxjj 11yj yj1(xxj1)(yj yj1)/(xj xj1)
缺点：I(x)连续，但不光滑，精度较低,仅在 hm 1jan{xhj xj xj1}足够小才能较好。的逼近
ss((xxn0
) )
f f
( x0 ) (xn )
m0 mn
数值分析三次样条插值
若取等距节点 hi = h, i = 1,…, n –1
i
h h
h
1 2
i
1 i
1 2
di
6 2h
yi 1
2 yi h
yi 1
3 h3
( yi1
2 yi
yi1 )
i 1, 2,, n
例1. 对于给定的节点及函数值
k 0123 xk 1 2 4 5 f (xk ) 1 3 4 2 求满足自然边界条件S(x0 ) S(xn ) 0的三次样条 插值函数S(x),并求f (3)的近似值
Mi1
( x xi )2 2hi 1
yi1 hi 1
yi
hi 1 6
( M i 1
Mi )
于是
Si( xi )
hi 3
Mi
yi
yi1 hi
hi 6
M i 1
Si1( xi )
hi 1 3
Mi
yi1 hi 1
yi
hi 1 6
M i 1
解： 由M关系式
k
hk
hk hk 1
k
hk 1 hk hk 1
1 k
1
2 3
1
1 3
2
1 3
2
2 3
di
6
yi1 hi1
yi
yi
yi hi
1
hi hi1 6 f [ xi1, xi , xi1]
三次样条插值算法详解
如果S(x)是f (x)的三次样条插值函数,则其必满足
插值条件： 连续性条件：
一阶导数连续条件：
二阶导数连续条件：
S(x j ) y j , j 0,1,, n
lim
xx j
S(x)
S(xj )
yj,
j
1,, n
1
lim
xx j
S ( x)
S(x j
)
mj
,
j
1,, n
1
lim
xx j
S
(
x)
S(
S(x)
(3x
3
16 x 2
27 x
14)
15
(x3 8x2 21x 18) 15
0 x 1 1 x 2
2 x3
10
三次样条插值函数的求法
通常有三转角法、三弯矩法、B样条基函数法。
这三种方法的基本思想是类似的，都是通过待定 某些参数来确定插值函数，但肯定不是待定4n个参
数。而是利用已知条件将待定参数减小到最少。
第一边界条件：由区间端点处的一阶导数给出即
s3 (x0 ) m0 f (x0 ), s3 (xn ) mn f (xn ),
6
第二边界条件：由区间端点处的二阶导数给出即
s3(x0 ) M 0 f (x0 ),
s3(
xn
)
Mn
f (xn ),
特殊情况为自然边界条件：
由区间端点处的二阶导数恒为0给出即
化为矩阵形式
17
2 1
2
2
2
m1 g1 1m0
m2
g2
3 2 3 4 2
m3
g3
n2 2 n2 mn2
数值分析课程设计--三次样条插值
《数值分析》课程设计三次样条插值算法院（系）名称信息工程学院专业班级 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）。
数值分析三次样条插值函数
数值分析三次样条插值函数【问题】对函数f x =ex, x∈[0,1]构造等距节点的三次样条插值函数，对以下两种类型的样条函数1. 三次自然样条2. 满足S′ 0 =1,S′ 1 =e的样条并计算如下误差：max{ f x1 −S x1 ,i=1,…,N} i−i−i这里xi−1为每个小区间的中点。
对N=10,20,40比较以上两组节点的结果。
讨论你的结果。
【三次样条插值】在每一个区间[t1,t2],…,[tn−1,tn]上，S都是不同的三次多项式，我们把在[ti−1,ti]上表示S的多项式记为Si，从而，S0 x x∈[t0,t1]∈[t1,t2] S x = S1 x x…Sn−1 x x∈[tn−1,tn]通过在节点处函数值、一阶导数和二阶导数的连续性可以得到：Si−1 ti = yi= Si ti 1≤i≤ n−1Si−1′ ti = Si′ tix→ti+limS′′ x =zi=limS′′(x) x→ti−再给定z0和zn 的值就构成了4n个条件，而三次样条插值函数共4n个系数，故可以通过这4n个条件求解三次样条函数的系数，从而求得该三次样条插值函数。
特别的，当z0=zn=0 时称为自然三次样条。
文本预览：一、自然三次样条插值【自然三次样条插值算法】1.由上面的分析可知，求解三次样条函数实际上就是求解一个矩阵：u 1h 1h1u2h2h2u3…v1 z1 v2 z2 z3=v3 … z…hn−2 n−2 vn−2 z vn−1 un−1 n−1ih3…hn−3un−2hn−26…其中hi=ti+1−ti，ui=2(hi+hi−1)，ui=h(yi+1−yi)，vi=bi−bi−1 所以自然三层次样条插值的算法就是在得到端点的函数值，一次导数值和二次导数值，然后根据上述求解矩阵得到v，代入自然三次样条的表达式即可。
2.根据题目中所给出的误差估计，计算在区间中点处的最大误差。
【实验】通过Mathematica编写程序得到如下结果：N=101. 计算得到zi的值为：由此可以得到各个区间的自然三次样条插值函数。
三次样条插值法
三次样条插值法根据李庆阳的《数值分析》这本教材中的讲解编写的程序，使用的是第一边界条件，用追赶法求解了M矩阵。
为了调用方便，我将整个函数所有的信息构造成一个结构体，输入插值点的坐标和数量，定义边界条件后，将这个结构体的指针作为参数传给Spline3（）函数，就完成了函数计算，计算结果也存储在该结构体中。
程序如下：/*=================================================================== ====*///=====================三次样条插值的函数S(x)实现=============================// 创建人：汪雅楠// 北京交通大学 QQ312818820// 说明：根据研究生教材《数值分析》（李庆杨）第51页~第56页编写/* 初始条件： 1. 已知两端的一阶导数值2. 已知两端的二阶导数值3. 周期样条函数### 此函数选择1条件 ###函数建立： 1. 设定样条函数S(x)的一阶导数为变量ki,用分段三次Hermitte差值2. 设定样条函数S(x)的二阶导数为变量Ki,用分段积分### 此函数选择2方法 ###矩阵求解：追赶法求解严格三对角占优矩阵{M},根据教材第195页编写*//*=================================================================== ====*/#include <stdio.h>///////////////////////////////////////////////////////////////////// ///////////#define MAXNUM 50 //定义样条数据区间个数最多为50个typedef struct SPLINE //定义样条结构体，用于存储一条样条的所有信息{ //初始化数据输入float x[MAXNUM+1]; //存储样条上的点的x坐标，最多51个点float y[MAXNUM+1]; //存储样条上的点的y坐标，最多51个点unsigned int point_num; //存储样条上的实际的点的个数float begin_k1; //开始点的一阶导数信息float end_k1; //终止点的一阶导数信息//float begin_k2; //开始点的二阶导数信息//float end_k2; //终止点的二阶导数信息//计算所得的样条函数S(x)float k1[MAXNUM+1]; //所有点的一阶导数信息float k2[MAXNUM+1]; //所有点的二阶导数信息//51个点之间有50个段，func[]存储每段的函数系数float a3[MAXNUM],a1[MAXNUM];float b3[MAXNUM],b1[MAXNUM];//分段函数的形式为Si(x) = a3[i] * {x(i+1) - x}^3 + a1[i] * {x(i+1) - x} +// b3[i] * {x - x(i)}^3 + b1[i] * {x - x(i)}//xi为x[i]的值，xi_1为x[i+1]的值}SPLINE,*pSPLINE;typedef int RESULT; //返回函数执行的结果状态，下面为具体的返回选项#ifndef TRUE#define TRUE 1#endif#ifndef FALSE#define FALSE -1#endif#ifndef NULL#define NULL 0#endif#ifndef ERR#define ERR -2#endif///////////////////////////////////////////////////////////////////// //////////////*=================================================================== ============*** 函数名称：Spline3()*** 功能说明：完成三次样条差值，其中使用追赶法求解M矩阵*** 入口参数：(pSPLINE)pLine 样条结构体指针pLine中的x[],y[],num,begin_k1,end_k1*** 出口参数：(pSPLINE)pLine 样条结构体指针pLine中的函数参数*** 返回参数：返回程序执行结果的状态TRUE or FALSE===================================================================== ===========*/RESULT Spline3(pSPLINE pLine){float H[MAXNUM] = {0}; //小区间的步长float Fi[MAXNUM] = {0}; //中间量float U[MAXNUM+1] = {0}; //中间量float A[MAXNUM+1] = {0}; //中间量float D[MAXNUM+1] = {0}; //中间量float M[MAXNUM+1] = {0}; //M矩阵float B[MAXNUM+1] = {0}; //追赶法中间量float Y[MAXNUM+1] = {0}; //追赶法中间变量int i = 0;////////////////////////////////////////计算中间参数if((pLine->point_num < 3) || (pLine->point_num > MAXNUM + 1)){return ERR; //输入数据点个数太少或太多}for(i = 0;i <= pLine->point_num - 2;i++){ //求H[i]H[i] = pLine->x[i+1] - pLine->x[i];Fi[i] = (pLine->y[i+1] - pLine->y[i]) / H[i]; //求F[x(i),x(i+1)] }for(i = 1;i <= pLine->point_num - 2;i++){ //求U[i]和A[i]和D[i]U[i] = H[i-1] / (H[i-1] + H[i]);A[i] = H[i] / (H[i-1] + H[i]);D[i] = 6 * (Fi[i] - Fi[i-1]) / (H[i-1] + H[i]);}//若边界条件为1号条件，则U[i] = 1;A[0] = 1;D[0] = 6 * (Fi[0] - pLine->begin_k1) / H[0];D[i] = 6 * (pLine->end_k1 - Fi[i-1]) / H[i-1];//若边界条件为2号条件，则//U[i] = 0;//A[0] = 0;//D[0] = 2 * begin_k2;//D[i] = 2 * end_k2;/////////////////////////////////////////追赶法求解M矩阵B[0] = A[0] / 2;for(i = 1;i <= pLine->point_num - 2;i++){B[i] = A[i] / (2 - U[i] * B[i-1]);}Y[0] = D[0] / 2;for(i = 1;i <= pLine->point_num - 1;i++){Y[i] = (D[i] - U[i] * Y[i-1]) / (2 - U[i] * B[i-1]);}M[pLine->point_num - 1] = Y[pLine->point_num - 1];for(i = pLine->point_num - 1;i > 0;i--){M[i-1] = Y[i-1] - B[i-1] * M[i];}//////////////////////////////////////////计算方程组最终结果for(i = 0;i <= pLine->point_num - 2;i++){pLine->a3[i] = M[i] / (6 * H[i]);pLine->a1[i] = (pLine->y[i] - M[i] * H[i] * H[i] / 6) / H[i];pLine->b3[i] = M[i+1] / (6 * H[i]);pLine->b1[i] = (pLine->y[i+1] - M[i+1] * H[i] * H[i] / 6) /H[i]; }return TRUE;}///////////////////////////////////////////////////////////////////// /////////////SPLINE line1;pSPLINE pLine1 = &line1;///////////////////////////////////////////////////////////////////// /////////////main(){line1.x[0] = 27.7;line1.x[1] = 28;line1.x[2] = 29;line1.x[3] = 30;line1.y[0] = 4.1;line1.y[1] = 4.3;line1.y[2] = 4.1;line1.y[3] = 3.0;line1.point_num = 4;line1.begin_k1 = 3.0;line1.end_k1 = -4.0;Spline3(pLine1);return 0;}///////////////////////////////////////////////////////////////////// /////////////。
三次样条插值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++代码实现。
数值分析(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 ];
数值分析
第五节 样条插值
样条是绘图员用于描绘光滑曲线的一种机 械器件,它是一些易弯曲材料制成的窄条或棒条. 在绘制需要通过某点的光滑曲线时,对它在这些 点的位置上“压铁”,它就被强制通过或接近图 表上确定的描绘点.“样条函数”这个术语意在 点出这种函数的图象与机械样条画出的曲线很 象.
数值分析实验报告-插值、三次样条
实验报告：牛顿差值多项式&三次样条... . (1)问题：在区间［-1,1］上分别取n=10、20用两组等距节点对龙格函数f (x)---作多项式插25 x 2值及三次样条插值对每个n值，分别画出插值函数矽(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 A10+494.91*x A8-9.5065e-14*x A7-381.43*x A6-8.504e-14*x A5+123.36*x A4+2.0202e-14*x A3-16.855*x A2-6.6594e-16*x+1.0并且这里也能得到该牛顿插值多项式的在［-1，1］上的图形，并和原函数进行对比（见Fig.1）。
数值分析2 6三次样条插值
4m h j 1
j
6 h2j1
(
yj
y j1 )
由条件
S( x j 0) S( x j 0) ( j 1,...,n 1)
可得
1 h j 1
m
j
1
2(
1 h j 1
1 hj
)m
j
1 hj
m
j1
3(
y
j
1 h2j
yj
y
j
y h2j1
j1
)
( j 1,...,n 1)
进一步简化为
jm j1 2m j jm j1 g j ( j 1,...n 1)
yk k (x) yk 1 k1(x)
其中
kk1((xx))(1(122xxkxxk1xxxkkxkk11)()(xxxkkx1xxkxkkx11k))22
k(
k 1
(
x) x)
(x (x
xk
)(
x xk
xk1 xk1
)2
xk
1
)(
x x xk1
k
x
k
)2
一、 三次样条的产生和背景
2.三次样条插值函数的定义
三次样条函数 +
S(xi) = yi
3.求解三次样条插值函数的已知条件数和 未知条件数
未知参数个数
4n
已知条件个数
插值条件：
n+1
S(x)∈C2[a,b] ：3(n-1)
共 计：
4n-2
缺少条件，通常在插值区间的端点给出，称 为边界条件。
4.常用的三种边界条件
1°已知两端的一阶导数值，即：
周期样条
S( x0 0) S( xn 0)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算方法作业实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。
实验函数:dt ex f xt ⎰∞--=2221)(π实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。
实验4.5 三次样条差值函数的收敛性实验目的：多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一算法描述：拉格朗日插值：错误!未找到引用源。
其中错误!未找到引用源。
是拉格朗日基函数，其表达式为：()∏≠=--=ni j j j i ji x x x x x l 0)()(牛顿插值：))...()(](,...,,[....))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N其中⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧--=--=--=-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i ji j i j i三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-1,x i ]上是三次多项式，其在此区间上的表达式如下：],[),6()6(]6)([6)(6)()(111113131i i ii i i i i i i i i i i i i i i i i i x x x h yM h M h h y x M M h h y y h x x Mi h x x M x S -------∈-+-+---+-+-=式中Mi=)(i x S ''.因此，只要确定了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 车门曲线。