Matlab Spline 三次样条插值多项式表达式问题
如何运用MATLAB 三次样条插值的问题，今天做作业，突然想用Matlab搞搞。
题目如下：
清华大学出版社的《数值分析（第5版）》
P49,20题。
x=[0.25 0.3 0.39 0.45 0.53];
y=[ 0.5 0.5477 0.6245 0.6708 0.7280 ]
pp=csape(x,y,'second',[0,0.0]);
disp(pp.coefs);
其中COEFS的含义是在Xi-Xi+1区间上的多项式是，例如COEFS数组第一行的意思是在X=0.25到X=0.3的区间上时表达式是-6.2652*（X-0.25)^3+0.9697*(X-0.25)^1+0.5;
-6.2652 0.0000 0.9697 0.5000
1.8813 -0.9398 0.9227 0.5477
-0.4600 -0.4318 0.7992 0.6245
2.1442 -0.5146 0.7424 0.6708
关于csape的用法引用自：https://www.360docs.net/doc/7c14636904.html,/ck436/blog/item/6fe40c46400d3c046b63e52b.html
csape，是计算在各种边界条件下的三次样条插值。
pp = csape(x,y,conds)
其中conds主要有以下的选项variational（自然边界条件，首末点二阶导数均为0），second （指定首末点的二阶导数），periodic（周期性边界条件，首末点的0~2阶导数相等），complete （给定导数情况，默认）
function pp = csape(x,y,conds,valconds)
%pp=csape(x,y,'变界类型','边界值'),生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量
%边界类型可为:'complete',给定边界一阶导数.
% 'not-a-knot',非扭结条件,不用给边界值.
% 'periodic',周期性边界条件,不用给边界值.
% 'second',给定边界二阶导数.
% 'variational',自然样条(边界二阶导数为0)
% .
%例考虑数据
% x | 1 2 4 5
% ---|-------------
% y | 1 3 4 2
%边界条件S''(1)=2.5,S''(5)=-3,
% x=[1 2 4 5];y=[1 3 4 2];
% pp=csape(x,y,'second',[2.5,-3]);pp.coefs % xi=1:0.1:5;yi=ppval(pp,xi);
% plot(x,y,'o',xi,yi);
三次样条插值---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] 画图验证：
运用matlab建立三次样条插值函数
（1）编写三条样条插值函数程序如下： x=[1 4 9 16 25 36 49 64 81]; y=[1 2 3 4 5 6 7 8 9]; n=length(x); lamda(1)=1; miu(n)=1; h=diff(x); df=diff(y)./diff(x); d(1)=6*(df(1)-1/2)/h(1); d(n)=6*(0.5*81^-0.5-df(n-1))/h(n-1); for j=2:n-1 lamda(j)=h(j)/(h(j-1)+h(j)); miu(j)=h(j-1)/(h(j-1)+h(j)); d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j)); end miu=miu(2:end); u=diag(miu,-1);r=diag(lamda,1);a=diag(2*ones(1,n)); A=u+r+a; %求出矩阵形式的线性方程组 M=inv(A)*d'; %求出M值 syms g for j=1:n-1 s(j)=M(j)*(x(j+1)-g)^3/(6*h(j))+M(j+1)*((g-x(j))^3/(6*h(j)))+(y(j)-M( j)*h(j)^2/6)*(x(j+1)-g)/h(j)+(y(j+1)-M(j+1)*h(j)^2/6)*(g-x(j))/h(j); end format rat for j=1:n-1 S(j,:)=sym2poly(s(j)); %三条样条插值函数 end %生成三次样条插值函数图象 for j=1:n-1 x1=x(j):0.01:x(j+1); y1=polyval(S(j,:),x1); plot(x1,y1,x,y,'o'); title('spline 三次样条插值函数图象'); xlabel('x'); ylabel('y'); grid on; hold on; end
MATLAB三次样条插值之三弯矩法
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;
三次样条插值方法的应用
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
三次样条插值的MATLAB实现
MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤
三次样条插值的Matlab实现(自然边界和第一边界条件)
（第一边界条件）源代码：function y=yt1(x0,y0,f_0,f_n,x)_____________(1) %第一类边界条件下三次样条插值； %xi所求点； %yi所求点函数值； %x已知插值点； %y已知插值点函数值； %f_0左端点一次导数值； %f_n右端点一次导数值； n = length(x0); z = length(y0); h = zeros(n-1,1); k=zeros(n-2,1); l=zeros(n-2,1); S=2*eye(n); fori=1:n-1 h(i)= x0(i+1)-x0(i); end fori=1:n-2 k(i)= h(i+1)/(h(i+1)+h(i)); l(i)= 1-k(i);
end %对于第一种边界条件： k = [1;k];_______________________(2) l = [l;1];_______________________(3) %构建系数矩阵S： fori = 1:n-1 S(i,i+1) = k(i); S(i+1,i) = l(i); end %建立均差表： F=zeros(n-1,2); fori = 1:n-1 F(i,1) = (y0(i+1)-y0(i))/(x0(i+1)-x0(i)); end D = zeros(n-2,1); fori = 1:n-2 F(i,2) = (F(i+1,1)-F(i,1))/(x0(i+2)-x0(i)); D(i,1) = 6 * F(i,2); end %构建函数D： d0 = 6*(F(1,2)-f_0)/h(1);___________(4)
三次样条插值作业题
例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)(π x 0.0 0.1 0.2 0.3 0.4 F(x) 0.5000 0.5398 0.5793 0.6179 0.7554 求f(0.13)和f(0.36)的近似值 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考 实验名称 实验 4.3三次样条插值函数（P126） 4.5三次样条插值函数的收敛性（P127） 实验时间 姓名 班级 学号 成绩
虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 段数据如下： k x 0 1 2 3 4 5 6 7 8 9 10 k y 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 k y ' 0.8 0.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 非常类似前面的三弯矩法，这里的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) % 三转角样条插值 给定数据表如下： 试求三次样条插值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 (第一边界条件源代码: function y=yt1(x0,y0,f_0,f_n,x _____________(1 %第一类边界条件下三次样条插值; %xi 所求点; %yi所求点函数值; %x 已知插值点; %y 已知插值点函数值; %f_0左端点一次导数值; %f_n右端点一次导数值; n = length(x0; z = length(y0; h = zeros(n-1,1; k=zeros(n-2,1; l=zeros(n-2,1; S=2*eye(n; fori=1:n-1 h(i= x0(i+1-x0(i; end fori=1:n-2 k(i= h(i+1/(h(i+1+h(i; l(i= 1-k(i; end %对于第一种边界条件: k = [1;k]; _______________________(2 l = [l;1]; _______________________(3 %构建系数矩阵 S : fori = 1:n-1 S(i,i+1 = k(i; S(i+1,i = l(i; end %建立均差表: F=zeros(n-1,2; fori = 1:n-1 F(i,1 = (y0(i+1-y0(i/(x0(i+1-x0(i; end D = zeros(n-2,1; fori = 1:n-2 F(i,2 = (F(i+1,1-F(i,1/(x0(i+2-x0(i; D(i,1 = 6 * F(i,2; end %构建函数 D : d0 = 6*(F(1,2-f_0/h(1; ___________(4 dn = 6*(f_n-F(n-1,2/h(n-1; ___________(5 D = [d0;D;dn]; ______________(6 m= S\D; %寻找 x 所在位置,并求出对应插值: fori = 1:length(x for j = 1:n-1 if (x(i<=x0(j+1&(x(i>=x0(j y(i =( m(j*(x0(j+1-x(i^3/(6*h(j+... (m(j+1*(x(i-x0(j^3/(6*h(j+... (y0(j-(m(j*h(j^2/6*(x0(j+1-x(i/h(j+... (y0(j+1-(m(j+1*h(j^2/6*(x(i-x0(j/h(j ; break; else continue; end end end (2 (自然边界条件源代码: 仅仅需要对上面部分标注的位置做如下修改 : 学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：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、第二类边界条件：给定函数在端点处的二阶导数，即 (){} 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); Ｍ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); (){} 21 ()(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 ，算法组织如下： 三次样条插值多项式 ——计算物理实验作业四 陈万物理学2013级 主程序： clear,clc; format rat x = [1,4,9,16,25,36,49,64]; y = [1,2,3,4,5,6,7,8]; f1 = ; fn = 1/16; [a,b,c,d,M,S] = spline(x,y,f1,fn); 子程序1： function [a,b,c,d,M,S]=spline(x,y,f1,fn) % 三次样条插值函数 % x是插值节点的横坐标 % y是插值节点的纵坐标 % u是插值点的横坐标 % f1是左端点的一阶导数 % fn是右端点的一阶导数 % a是三对角矩阵对角线下边一行 % b是三对角矩阵对角线 % c是三对角矩阵对角线上边一行 % S是插值点的纵坐标 n = length(x); h = zeros(1,n-1); deltay = zeros(1,n); miu = zeros(1,n-1); lamda = zeros(1,n-1); d = zeros(1,n-1); for j = 1:n-1 h(j) = x(j+1)-x(j); deltay(j) = y(j+1)-y(j); end % 得到h矩阵 for j = 2:n-1 sumh = h(j-1) + h(j); miu(j) = h(j-1) / sumh; lamda(j) = h(j) / sumh; d(j) = 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)))/sumh; end % 根据第一类边界条件，作如下规定 lamda(1) = 1; d(1) = 6*(deltay(1)/h(1)-f1)/h(1); miu(1) = 1; d(n) = 6*(fn-deltay(n-1)/h(n-1))/h(n-1); 告： Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数p chip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x); 2.>>xx = 0:.25:10; yy = interp1(x,y,xx); 3.>>plot(x,y,'kd',xx,yy) 复制代码 例2 1.>> year = 1900:10:2010; 2.>> product = [75.995 91.972 105.711 12 3.203 131.669 150.697 179.323 203.212 226.505 数值分析上机报告 题目：插值法 学号：201014924 姓名：靳会有 一、调用MATLAB内带函数插值 1、MATLAB内带插值函数列举如下： 2、取其中的一维数据内插函数（）为例，程序如下：其调用格式为： yi=interp1(x, y, xi) yi=interp1(x, y, xi, method) 举例如下： x=0:10:100 y=[40 44 46 52 65 76 80 82 88 92 110]; xi=0:1:100 yi=interp1(x,y,xi,'spline') 3、其他内带函数调用格式为： Interpft函数： y=interpft(x,n) y=interpft(x,n,dim) interp2函数： ZI=interp2(X, Y, Z, XI, YI)，ZI=imerp2(Z, ntimes) ZI=interp2(Z, XI, YI) ，ZI=interp2(X, Y, Z, XI, YI, method) interp3函数： VI=interp3(X,Y,Z,V,XI,YI,ZI) VI=interp3(V, ntimes) VI=interp3(V,XI,YI,ZI) VI=interp3(…, method) Interpn函数： VI=interpn(X1, X2, X3, …, V, Y1, Y2, Y3, …) VI=interpn(V, ntimes) VI=interpn(V, Yl, Y2, Y3, …) VI=interpn(…, method) Spline函数： yi=spline(x,y,xi) pp=spline(x,y) meshgrid函数： [X,Y]=meshgrid(x,y) [X,Y]=meshgrid(x) [X,Y,Z]=meshgrid(x,y,z) Ndgrid函数： [X1, X2, X3, …]=ndgrid(x1, x2, x3, …) [X1, X2, X3, …]=ndgrid(x) Griddata函数： ZI=griddata(x, y, z, XI, YI) [XI, YI, ZI]=griddata(x, y, z, xi, yi) […]=griddata(… method) 二、自编函数插值 1、拉格朗日插值法： 建立M 文件： function f = Language(x,y,x0) syms t l; if(length(x) == length(y)) n = length(x); else disp('x和y的维数不相等！'); return; %检错MATLAB三次样条插值之三转角法
试求三次样条插值S(X)
三次样条插值的Matlab实现(自然边界和第一边界条件)(精)
关于三次样条插值函数的学习报告(研究生)资料
matlab_牛顿插值法_三次样条插值法
三次样条插值自然边界条件
MATLAB三次样条插值之三弯矩法
matlab 牛顿插值法 三次样条插值法
三次样条插值多项式matlab
Matlab中插值函数汇总和使用说明
MATLAB实现拉格朗日插值精编版