三次样条插值课程设计含代码和图
一、背景
L-插值（牛顿插值）
Hermite 插值 高次插值出现龙格现象 分段插值 分段线性插值在节点处不一定光滑
分段Hermite 插值 但导数值不容易提取（找到）
三次样条插值（先由函数值确定导数值，再由分段
Hermite 插值解决问题）
定义：
满足条件
如果函数：上给出一个划分
，在区间，
上的二次连续可微函数是区间设函数)(...][],[)(110x s b
x x x x a b a b a x f n n =<<<<=?- 次的多项式；
是一个不超过
上在每个小区间m x s n k x x k k )(),...,2,1](,[)1(1=+
阶的连续倒数。处具有（在节点1m ),...,2,1)3(-=n k x k );,...2,1,0()()(1n j x f x s j j ==）（
一、上机实验的问题和要求（需求分析）：
掌握三次样条插值法的上机编程运算。
掌握分析插值结果优劣的方法
实验内容
已知 x = [1 2 4 5 6 7];
y = [1 3 4 2 3 4];
试计算f （3）近似值，计算与精确值之间的误差，验证实际结果和理论 结论的一致性。
二、程序设计的基本思想，原理和算法描述：
三、主要程序代码或命令：
clear; clc; close all;
x = [1 2 4 5 6 7];
y = [1 3 4 2 3 4];
n = length(x);
A = zeros(n,n);
d = zeros(n,1);
m=zeros(n,1);
mu = zeros(n-2,1);
lambda = zeros(n-2,1);
h = zeros(n-1,1);
lambda = zeros(n-2,1);
h = zeros(n-1,1);
h(1) = x(2) - x(1);
for j = 2:n-1
h(j) = x(j+1) - x(j);
lambda(j-1) = h(j)/(h(j) + h(j-1));
mu(j-1) = 1 - lambda(j-1);
d(j) = 3*(mu(j-1)*(y(j+1)-y(j))/h(j)+lambda(j-1)*(y(j)-y(j-1))/h(j-1)); end
for j = 2:n-1
A(j,j) = 2;
A(j,j+1) = mu(j-1);
A(j+1,j) = lambda(j-1);
end
A(n,n) = 2;
A(1,1)=2;A(1,2)=1;A(n-1,n)=1;A(n,n)=2;
d(1)=3*(y(2)-y(1)/h(1));
d(n)=3*(y(n)-y(n-1))/h(n-1);
m = A\d;
C = zeros(n-1,3);
for j = 1:n-1
C(j,1) = m(j);
C(j,2) = (3*(y(j+1)-y(j))/h(j)-2*m(j)-m(j+1))/h(j);
C(j,3) = (m(j+1)+m(j)-2*(y(j+1)-y(j))/h(j))/(h(j)*h(j));
end
h_test = 0.01;
m = length(x_test);
y_test = zeros(m,1);
j = 1;
for i = 1:m
if x_test(i) <= x(j+1)
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
else
j = j + 1;
y_test(i) = y(j) + C(j,1)*(x_test(i) - x(j)) + C(j,2)*(x_test(i) - x(j))^2 + C(j,3)*(x_test(i) - x(j))^3;
end
plot(x,y,'*',x_test,y_test,'r')
legend('2??μ?úμ?','èy′??ùì?',2)
axis auto
x_nu = 3;
for j = 1:n-1
if x_nu > x(j) && x_nu <= x(j+1)
break;
else
continue;
end
end
y_nu = y(j) + C(j,1)*(x_nu - x(j)) + C(j,2)*(x_nu - x(j))^2 + C(j,3)*(x_nu - x(j))^3;
y_nu
四、调试和运行程序过程中产生的问题及采取的措施：
下标不统一，matlab从1开始
五、运行输出结果及分析：
三次样条插值代码
2 三次样条插值程序 三次样条插值利用方案二（求解固支样条或压紧样条） 按照要求要起点和终点的一阶导数值已知， 可得关于01,,.....,n M M M 的严格对角占优势的三对角方程组 然后利用三对角法（追赶法）解此线性方程组。 （1）编写M 文件，并保存文件名scfit.m % x,y 分别为n 个节点的横坐标和纵坐标值组成的向量 % dx0和dxn 分别为S 的导数在x0和xn 处的值，即m 0和m n n=length(x)-1; h=diff(x); d=diff(y)./h; a=h(2:n-1); b=2*(h(1:n-1)+h(2:n)); c=h(2:n); u=6*diff(d); b(1)=b(1)-h(1)/2; u(1)=u(1)-3*(d(1)-dx0); b(n-1)=b(n-1)-h(n)/2; u(n-1)=u(n-1)-3*(dxn-d(n)); %追赶法部分 for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1); end m(n)=u(n-1)/b(n-1); for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2))/b(k); end %求S K1，S K2，S K3，S K4 m(1)=3*(d(1)-dx0)/h(1)-m(2)/2; m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2; for k=0:n-1 00 ()S x m '=()n n S x m '=0011111111212212n n n n n n M d M d M d M d μλμλ----??????????????????????=??????????????????????????
三次样条插值---matlab实现
计算方法实验—三次样条插值 机电学院075094-19 苏建加 20091002764 题目：求压紧三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且一阶导 数边界条件S'(-3)=-1;S'（4）=1。 解：首先计算下面的值： 记 1--=j j j x x h ； 1++=j j j j h h h u ；1=+j j u λ ； ?? ????????---+=-++++-j j j j j j j j j j j h y y h y y h h x x x f 1111 111],,[ ；M j =)(''j x s ;],,[611+-=j j j j x x x f d ; h1=-2-(-3)=1;h2=1-(-2)=3;h3=4-1=3; u1=1/4;u2=3/6; d1=6/4*(3/3-(-2)/1)=4.5;d2=6/6*(-2/3-3/3)=-5/3; 由于边界条件S'(-3)=-1;S'（4）=1,得到如下 式子： d0=6/1*(-2/1-(-1))=-6； d3=6/3*(1-(-2)/3)=10/3; 所以得到4个含参数m0~m3 的线性代数方程组为： 2.0000 1.0000 0 0 m0 0.2500 2.0000 0.7500 0 m1 0 0.5000 2.0000 0.5000 m2 0 0 1.0000 2.0000 m3 利用matlab 求解方程得： m = -4.9032 3.8065 -2.5161 2.9247 所以 S1(x)=-0.8172*(-2-x)^3+ 0.6344*(x+3)^3+2.8172*(-2-x)-0.6344*(x+3) x ∈[-3，-2] S2(x)=0.2115*(1-x)^3 -0.1398*(x+2)^3- 1.9032*(1-x)+ 2.2581*(x+2) x ∈[-2，1] S3(x)=-0.1398*(4-x)^3+0.1625(x-1)^3+ 2.2581*(4-x)-1.1290*(x-1) x ∈[1，4] 化简后得：S1(x)=1.4516*x^3 + 10.6128*x^2 + 23.4836*x + 16.1288 x ∈[-3，-2] S2(x)=-0.3513x^3-0.2043x^2+1.8492x+1.7061 x ∈[-2，1] S3(x)=0.3023x^3-2.1651x^2+3.8108x+1.0517 x ∈[1，4] 画图验证：
三次样条插值函数
沈阳航空航天大学 数学软件课程设计 （设计程序） 题目三次样条插值函数 班级 / 学号 学生姓名 指导教师
沈阳航空航天大学 课程设计任务书 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学 班级学号姓名 课程设计题目三次样条插值函数 课程设计时间: 2010 年12月20日至2010 年12月31日 课程设计的内容及要求： 1.三次样条插值函数 给出函数在互异点处的值分别为。 (1)掌握求三次样条插值函数的基本原理； (2)编写程序求在第一边界条件下函数的三次样条插值函数； (3)在区间上取n=10，20，分别用等距节点对函数 作三次样条插值函数，利用(1)的结果画出插值函数的图形，并在该图形界面中同时画出的图形。 [要求] 1.学习态度要认真，要积极参与课程设计，锻炼独立思考能力； 2.严格遵守上机时间安排； 3.按照MATLAB编程训练的任务要求来编写程序； 4.根据任务书来完成课程设计论文； 5.报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”； 6.报告上交时间：课程设计结束时上交报告；
7.严谨抄袭行为。 指导教师年月日负责教师年月日学生签字年月日
沈阳航空航天大学 课程设计成绩评定单 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学课程设计题目三次样条插值函数 学号姓名 指导教师评语： 课程设计成绩 指导教师签字 年月日
目录 一正文 (1) 1问题分析 (1) 1.1 题目 (1) 1.2 分析 (1) 2 研究方法原理 (1) 2.1 求三次样条插值多项式，算法组织 (1) 3 算例结果 (3) 二总结 (7) 参考文献 (8) 附录 (9) 源程序: (9) 程序1 (9) 程序2 (10) 程序3 (12) 程序 4 (12)
三次样条插值、拉格朗日插值、herminte插值
三次样条插值： function s=spline(x0,y0,y2l,y2n,x) n=length(x0); km=length(x); a(1)=-0.5; b(1)=3*(y0(2)-y0(1))/(2*(x0(2)-x0(1))); for j=1:n-1 h(j)=x0(j+1)-x0(j); end for j=2:n-1 alpha(j)=h(j-1)/(h(j-1)+h(j)); beta(j)=3*((1-alpha(j))*y0(j)-y(j-1)/h(j-1)+alpha(j)*(y0(j+1)-y0(j))/h(j)); a(j)=-alpha(j)/(2+(1-alpha(j))*a(j-1)); b(j)=(beta(j)-(1-alpha(j))*b(j-1))/(2+(1-alpha(j))*a(j-1)); end m(n)=(3*(y0(n)-y0(n-1))/h(n-1)+y2n*h(n-1)/2-b(n-1))/(2+a(n-1)); for j=(n-1):-1:1 m(j)=a(j)*m(j+1)+b(j); end for k=1:km for j=1:(n-1) if ((x(k)>x0(j))&(x(k) MATLAB三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（1，2，…，n）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数 sanwanj得到方程的系数矩阵，即追赶法方程的四个向量参数，接下来调用 追赶法（在intersanwj函数中），得到三次样条分段函数系数因子，然后进行多项式合并得 到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改 点的值。附：追赶法程序 chase %%%%%%%%%%%%%% function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b) % 三弯矩样条插值 % 将插值点分两次输入，x0 y0 单独输入 % 边值条件a的二阶导数 y1a 和b的二阶导数 y1b，这里建议将y1a和y1b换成y2a和 y2b，以便于和三转角代码相区别 n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1); w=2*ones(n+1); h0=x(1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1)); % for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k)); end newv=[v;1]; newu=[1;u]; d0=6*((y(1)-y0)/h0-y1a)/h0; 例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表： 且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s 本算法求解出的三次样条插值函数将写成三弯矩方程的形式： ) ()6()() 6()(6)(6)(211123 13 1j j j j j j j j j j j j j j j j x x h h M y x x h h M y x x h M x x h M x s -- + -- + -+ -= +++++其中，方程中的系数 j j h M 6， j j h M 61+，j j j j h h M y )6(2- ， j j j 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) - 1 h(i) = IndVar(i + 1) - IndVar(i); end % 为向量μ赋值 数值计算方法作业 实验4.3 三次样条差值函数 实验目的： 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考 虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 算法描述： 拉格朗日插值： 错误!未找到引用源。 其中错误!未找到引用源。是拉格朗日基函数，其表达式为：() ∏ ≠=--=n i j j j i j i 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 j i j i j i 三样条插值： 所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a 实验报告：牛顿差值多项式&三次样条 问题：在区间[-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 all clc x1=-1:0.2:1; y1=1./(1+25.*x1.^2); n=length(x1); f=y1(:); for j=2:n for i=n:-1:j f(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1)); end end syms F x p ; F(1)=1;p(1)=y1(1); for i=2:n F(i)=F(i-1)*(x-x1(i-1)); p(i)=f(i)*F(i); end syms P P=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 on xlabel('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.0 202e-14*x^3-16.855*x^2-6.6594e-16*x+1.0 并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。 Fig.1 牛顿插值多项式（n=10）函数和原函数图形 从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时，这一现象将更加明显。 1.2 当n=20时： 对n=10的代码进行修改就可以得到n=20时的代码。将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x^20 - 1.0121e6*x^18 + 2.6193e-12*x^17 + 1.6392e6*x^16 + 2.248e-11*x^15 - 1.4429e6*x^14 - 4.6331e-11*x^13 + 757299.0*x^12 + 1.7687e-11*x^11 - 245255.0*x^10 + 2.1019e-11*x^9 + 49318.0*x^8 + 3.5903e-12*x^7 - 6119.2*x^6 - 1.5935e-12*x^5 + 470.85*x^4 + 1.3597e-14*x^3 - 24.143*x^2 - 1.738e-14*x + 1.0 同样的，这里得到了该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.2）。 摘要 本文细致的讲解了三次样条插值函数的产生及在实际中解决的问题，通过MATLAB的程序编写，可以将复杂的计算省去，直接的给出了三次样条插值的结果与实际结果间的误差，验证实际结果和理论值的一致性。避免了求解方程中的不必要计算，使求解效率得到显著的提高。 关键词插值函数三次样条插值 MATLAB 1 三次样条插值函数概论 当插值节点很多时,插值多项式的次数就会很高,这不仅增大了计算量,还会影响结果的精确度.虽然可以采用上述分段插值,但是主要缺点就是个分段接头处不光滑,插值函数的导数不连续,因此想构造这样的插值:既能分段的低次插值,又能保证接头处的光滑,就产生了三次样条插值函数. 1.1定义 设函数()f x 市区间[a,b]上的二次连续可微函数，在区间[a,b]上给处一个划分。设函数()f x 是区间[a,b]上的一个划分 011...n n a x x x x b -?=<<<<= 如果函数()S x 满足条件 （1）在每个小区间1[,]k k x x +（k=1,2,….,n ）上()S x 是一个部超过m 次的多项式。 （2）在节点k x （k=1,2,….,n ）处具有m-1阶的连续导数。 （3）()()(0,1,2,...) j j s x f x j n == 1.2三次样条差值函数的构造 由于三次样条插值我、函数s(x)的插值节点处的二阶导数存在，因此令各节点处的二阶导数为 ' ()(0,1,...,)k s x m k n == （1.01） 根据样条插值函数的定义，三次样条插值函数是s(x)在每一个小区间)1....,1,0](,[]1-=+n k x x k k 上市不超过三次的多项式。在每一个小区间 )1....,1,0](,[]1-=+n k x x k k 上，其二阶导数为线性函数，即 '' 11 11()k k k k k k k k x x x x s x m m x x x x ++++--=+-- （1.02） 对式（1.02）积分两次，则得到 k k k k k k k k k b x x a h x x m h x x m x s +-+-++=++)(6)(6)()(3 1 3 1 （1.03） 实验二插值法 1、实验目的： 1、掌握直接利用拉格郎日插值多项式计算函数在已知点的函数值；观察拉格郎日插值的龙格现象。 2、了解Hermite插值法、三次样条插值法原理，结合计算公式，确定函数值。 2、实验要求: 1)认真分析题目的条件和要求，复习相关的理论知识，选择适当的解决方案和算法； 2)编写上机实验程序，作好上机前的准备工作； 3)上机调试程序，并试算各种方案，记录计算的结果（包括必要的中间结果）； 4)分析和解释计算结果； 5)按照要求书写实验报告； 3、实验内容： 1) 用拉格郎日插值公式确定函数值；对函数f(x)进行拉格郎日插值，并对f(x)与插值多项式的曲线作比较。 已知函数表：（0.56160,0.82741）、（0.56280,0.82659）、（0.56401,0.82577）、（0.56521,0.82495）用三次拉格朗日插值多项式求x=0.5635时函数近似值。 2) 求满足插值条件的插值多项式及余项 1) 4、题目：插值法 5、原理： 拉格郎日插值原理： n次拉格朗日插值多项式为：L n (x)=y l (x)+y 1 l 1 (x)+y 2 l 2 (x)+…+y n l n (x) n=1时，称为线性插值， L 1(x)=y (x-x 1 )/(x -x 1 )+y 1 (x-x )/(x 1 -x )=y +(y 1 -x )(x-x )/(x 1 -x ) n=2时，称为二次插值或抛物线插值， L 2(x)=y (x-x 1 )(x-x 2 )/(x -x 1 )/(x -x 2 )+y 1 (x-x )(x-x 2 )/(x 1 -x )/(x 1 -x 2 )+y 2 (x -x 0)(x-x 1 )/(x 2 -x )/(x 2 -x 1 ) n=i时， Li= （X-X0）……(X-X i-1)(x-x i+1) ……(x-x n) （X-X0）……(X-X i-1)(x-x i+1) ……(x-x n) 6、设计思想： 拉格朗日插值法是根据n + 1个点x0, x1, ... x n(x0 < x1 < ... x n)的函数值f (x0), f (x1) , ... , f (x n)推出n次多項式p(x)，然后n次多項式p (x)求出任意的点x对应的函数值f (x)的算法。 7、对应程序： 1 ) 三次拉格朗日插值多项式求x=0.5635时函数近似值 #include"stdio.h" #define n 5 void main() { int i,j; float x[n],y[n]; float x1; float a=1; float b=1; float lx=0; printf("\n请输入想要求解的X：\n x="); scanf("%f",&x1); printf("请输入所有点的横纵坐标：\n"); for(i=1;i 给定数据表如下： 试求三次样条插值S（X），并满足条件： i)S’(0.25)=1.0000, S’(0.53)-0.6868; ii) S”(0.25)= S”（0.53）=0； 解： 由给定数据知： h0 =0.3-0.25 - 0.05 , h 1=0.39-0.30-0.09 h 2=0.45-0.39-0.06, h 3=0.53-0.45-0.08 由μ i=h i/(h i1+h i), λ i= h i/(h i1+h i) 得: μ1= 5/14 ; λ 1= 9/14 μ2= 3/5 ; λ 2= 2/5 μ3= 3/7 ; λ 3=4/7 0.25 0.5000 ﹨ ﹨ 1.0000 ∕﹨ 0.25 0.5000 ∕ -0.9200-f[x 0,x 0, x 1 ] ﹨∕ 0.9540 ∕﹨ 0.30 0.5477 -0.7193-f[x 0,x 1,x 2 ] ﹨∕ 0.8533 ∕﹨ 0.39 0.6245 -0.5440-f[x1,x2,x 3 ] ﹨∕ 0.7717 ∕﹨ 0.45 0.6708 -0.4050-f[x 2,x 3,x 4 ] ﹨∕ 0.7150 ∕﹨ 0.53 0.7280 -0.3525-f[x 3,x 4,x 5 ] ﹨∕ 0.6868 ∕ 0.53 0.7280 i)已知一节导数边界条件，弯矩方程组 ┌┐┌┐ │ 2 1 │┌M 0 ┐│-0.9200 ︳ ︳5/14 2 9/14 ︳︳M ︳︳-0.7193 ︳ 1 ︳3/5 2 2/5 ︳︳M 2 ︳_6 ︳-0.5440︳ ︳ 3/7 2 4/7 ︳︳M ︳︳-0.4050 ︳ 3 MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤ 学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：152111033 姓名：刘楠楠 样条函数： 由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。 一、三次样条函数的定义： 对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<< <<=，若 函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在 [,]a b 上的三次样条函数。 二、三次样条函数的确定： 由定义可设：101212 1(),[,] (),[,]()(),[,] n n n s x x x x s x x x x s x s x x x x -∈??∈?=???∈?其中()k s x 为1[,]k k x x -上的三次 多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n = 由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++= ''''1()(),(1 ,2,,1)k k k k s x s x k n -+ +==-， 已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。 1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f == 2、第二类边界条件：给定函数在端点处的二阶导数，即 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 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。 Matlab 代码如下： function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second (){} 2 1 ()(11),5,10,20: 1252 1()1,(0,1,2,,)()2,(0,1,2,,)() ()2 35,20:1100 (i i i i n n k k k Newton f x x n x f x x i i n f x n x y i n Newton N x S x n x k y f x = -≤≤=+=-+====-+ = 题目：插值多项式和三次样条插值多项式。已知对作、计算函数在点处的值；、求插值数据点 的插值多项式和三次样条插值多项式；、对计算和相应的函数值),()() (1,2,,99)4:()max ()()max ()n k n k n k n k n k n k k k N x S x k E N y N x E S y S x ==-=- 和； 、计算，； 解释你所得到的结果。 算法组织: 本题在算法上需要解决的问题主要是：求出第二问中的Newton 插值多项式 )(x N n 和三次样条插值多项式()n S x 。如此，则第三、四问则迎刃而解。计算两 种插值多项式的算法如下： 一、求Newton 插值多项式)(x N n ，算法组织如下： Newton 插值多项式的表达式如下： )())(()()(110010--???--+???+-+=n n n x x x x x x c x x c c x N 其中每一项的系数c i 的表达式如下： 1102110) ,,,(),,,(),,,(x x x x x f x x x f x x x f c i i i i i -???-???= ???=- 根据i c 以上公式，计算的步骤如下： ?? ??? ?? ?????+??????? ???????????----) ,,,,(1) ,,,(),,,,(),(,),,(2)(,),(),(11101111011010n n n n n n n n x x x x f n x x x f x x x f n x x f x x f x f x f x f 、计算、计算、计算、计算 二、求三次样条插值多项式)(x S n ，算法组织如下： 例：已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 0.25 0.3 0.39 0.45 0.53 0.5000 0.5477 0.6245 0.6708 0.7280 求解，其中边界条件为. 1）三次样条插值自然边界条件源程序： function s=spline3(x,y,dy1,dyn) %x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导 m=length(x);n=length(y); if m~=n error('x or y输入有误') return end h=zeros(1,n-1); h(n-1)=x(n)-x(n-1); for k=1:n-2 h(k)=x(k+1)-x(k); v(k)=h(k+1)/(h(k+1)+h(k)); u(k)=1-v(k); end g(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1; g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn; for i=2:n-1 g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); end for i=2:n-1; A(i,i-1)=v(i-1); A(i,i+1)=u(i-1); end A(n,n-1)=1; A(1,2)=1; A=A+2*eye(n); M=zhuigf(A,g); %调用函数，追赶法求M fprintf('三次样条（三对角）插值的函数表达式\n'); syms X; for k=1:n-1 fprintf('S%d--%d:\n',k,k+1); s(k)=(h(k)+2*(X-x(k)))./h(k).^3.*(X-x(k+1)).^2.*y(k)... +(h(k)-2*(X-x(k+1)))./h(k).^3.*(X-x(k)).^2.*y(k+1)... +(X-x(k)).*(X-x(k+1)).^2./h(k).^2*M(k)+(X-x(k+1)).*... (X-x(k)).^2./h(k).^2*M(k+1); end s=s.'; s=vpa(s,4); %画三次样条插值函数图像 for i=1:n-1 X=x(i):0.01:x(i+1); st=(h(i)+2*(X-x(i)))./(h(i)^3).*(X-x(i+1)).^2.*y(i)... +(h(i)-2.*(X-x(i+1)))./(h(i)^3).*(X-x(i)).^2.*y(i+1)... +(X-x(i)).*(X-x(i+1)).^2./h(i)^2*M(i)+(X-x(i+1)).*... (X-x(i)).^2./h(i)^2*M(i+1); plot(x,y,'o',X,st); hold on End plot(x,y); grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %调用的函数： %追赶法 function M=zhuigf(A,g) n=length(A); L=eye(n); U=zeros(n); for i=1:n-1 U(i,i+1)=A(i,i+1); end U(1,1)=A(1,1); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i); end Y(1)=g(1); for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1); end M(n)=Y(n)/U(n,n); for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i); 非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。 %%%%%%%%%%%%%%%%%%% function [newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b) % 三转角样条插值 % 将插值点分两次输入，x0 y0 单独输入 % 边值条件a的一阶导数 y1a 和b的一阶导数 y1b n=length(x);m=length(y); if m~=n error('x or y 输入有误，再来'); end v=ones(n-1,1); u=ones(n-1,1); d=zeros(n-1,1); w=2*ones(n-1,1); h0=x(1)-x0; h=zeros(n-1,1); for k=1:n-1 h(k)=x(k+1)-x(k); end v(1)=h0/(h0+h(1)); u(1)=1-v(1); d(1)=3*(v(1)*(y(2)-y(1))/h(1)+u(1)*((y(1)-y0))/h0); % for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k)); u(k)=1-v(k); d(k)=3*(v(k)*(y(k+1)-y(k))/h(k)+u(k)*(y(k)-y(k-1))/h(k-1)); end d(1)=d(1)-u(1)*y1a; d(n-1)=d(n-1)-v(n-1)*y1b; newv=v(1:n-2,:); newu=u(2:n-1,:); %%%%%%%%%%%% function intersanzhj(x,y,x0,y0,y1a,y1b) % 三转角样条插值 ＭATLＡＢ三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（１，２，…,ｎ）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数saｎwanｊ得到方程的系数矩阵，即追赶法方程的四个向量参数,接下来调 用追赶法(在inteｒｓaｎwj函数中），得到三次样条分段函数系数因子,然后进行多项式合并 得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算 改点的值。附:追赶法程序cｈase ％%%％％%%%%%%%％% fｕncｔiｏｎ [nｅwv,w,newu，neｗd]=sanwj（ｘ,y,x０,y0,y1ａ,y1b）?％三弯矩样 条插值?％将插值点分两次输入，x0y０单独输入?% 边值条件a的二阶导数 y1a 和ｂ 的二阶导数ｙ1b，这里建议将y１ａ和y1b换成y2a和ｙ2b,以便于和三转角代码相区别 ?n=length(x);m=leｎgth(y); ｉf m~=n?errｏｒ(＇x or y 输入有误,再来'）; end?v=ｏneｓ(ｎ-１，1);u=ｏnes(n－１,1）;d=zeros（n－１,１）;?w=２*ｏ nes(n+１)；?h0=x（１）-x0;?h=zerｏs(n-１,1); for k＝1:n-１?ｈ(k)＝x（k+1)－x（k);?enｄ ｖ(1)＝h0/(ｈ0+h(1）); u(1)＝1－v（1）; ｄ（１)＝6*((y（2）-y(1）)/h（1)-(y(1)-y０)/ｈ0）/(h0+ｈ(１）);?% ｆoｒ k＝2:ｎ-1?ｖ(ｋ）＝h(ｋ-１)/(ｈ(k-１)+h(ｋ)）;?u(k)=1－v(k);?d（k)＝ 6*((y(ｋ+１)-ｙ（ｋ)）/h（k）-(ｙ(k）-y(ｋ-1）)/ｈ(k－1))／(h（k-1)+h(k)); ｅnd newv=[v;1］;?newu=［1;u]； d0=6＊（(y(1)-y0)/h0-y1a)/h0; d(n)=6＊（y1ｂ-(ｙ(n）－ｙ(n-1))／ｈ(n-1))／h(ｎ-1）; newd=［d0;d]; %%%%%%％％%％％％ ｆunctｉon iｎｔersａnｗj（x，y,ｘ0,y0,y1a，y1b) %三弯矩样条插值?%第一部分?n=leｎgtｈ(ｘ）;ｍ＝length(y); iｆ m~＝n?eｒror（'ｘｏｒy 输入有误，再来'); enｄ?%重新定义h?h=zｅｒｏs(n,1）; h(1)=x（1）-ｘ０; fｏr k=2:n h（k)=ｘ(ｋ)－x(k-1)；?end ％sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x，ｙ,ｘ0,ｙ0,y１ａ,y1b);MATLAB三次样条插值之三弯矩法
三次样条插值作业题
数值分析作业-三次样条插值
数值分析实验报告-插值、三次样条Word版
计算方法三次样条插值课程设计
插值法实验报告
试求三次样条插值S(X)
(精选)三次样条插值的MATLAB实现
关于三次样条插值函数的学习报告(研究生)资料
三次样条插值方法的应用
matlab_牛顿插值法_三次样条插值法
三次样条插值自然边界条件
MATLAB三次样条插值之三转角法
MATLAB三次样条插值之三弯矩法