MATLAB三次样条插值之三弯矩法
数值分析课程设计报告书三次样条插值的三弯矩法
数值分析课程设计报告书院系名称：学生姓名：专业名称：班级：时间：实验一 三次样条插值的三弯矩法一、实验目的已知数据i x ,()i i y f x =,0,,i n =及边界条件()n j x y j j 1,0),(2=,求)(x f 的三次样条插值函数)(x S .要求输出用追赶法解出的弯矩向量0[,,]n M M M =及()(),0,,,0,1,2k i S t i m k ==的值.画出)(x S y =的图形,图形中描出插值点(,)i i x y 及(,())i i t S t 分别用‘o ’和‘*’标记.二、实验原理1.用追赶法求解第二类边界条件的三弯矩方程:0010012111121111[,,]21[,,]26[,,]212[,,]n n n n n n n n n n f x x x M f x x x M M f x x x M f x x x μλμλ------⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 其中1111,,j jj j j j j j j j j h h h x x h h h h μλ-+--===-++.2.得出样条函数表达式:332211111()()()()()6666j j j j j j j j j j j j j j j jx x x x M h x x M h x x S x M M y y h h h h +++++----=++-+-. 3.计算(k)(),0,,,0,1,2i S t i m k ==.三、实验结果所用数据:x=[-2.223,-1.987,-1.8465,-1.292,-1.2266,-1.1056,-0.8662,-0.6594,-0.2671,-0.0452,0.5385,1.2564,1.4398,1.5415,1.7646,1.9678,2.236];y=[0.83995,1.1696,1.3141,1.6992,1.7312,1.7847,1.8708,1.9262,1.9881,1.9997,1.9511,1.7169,1.618,1.5543,1.3871,1.191,0.81662];d2s1= -4.5000;d2sn= -4.8967; %第二种边界条件t=[-2.223,-1.9443,-1.6656,-1.3869,-1.1083,-0.82956,-0.55088,-0.27219,0.0065,0.28519,0.56387,0.84256,1.1212,1.3999,1.6786,2.236]; ;(指定计算点)计算结果:-2.5-2-1.5-1-0.500.51 1.52 2.50.811.21.41.61.82四、实验分析通过实验结果我们，知道三弯矩法求出满足初始条件的三次样条函数，与其他插值函数的构造相比，三次样条插值法的计算量要小得多。
matlab三次样条插值例题解析
文章标题：深度解析Matlab三次样条插值1. 前言在数学和工程领域中，插值是一种常见的数值分析技术，它可以用来估计不连续数据点之间的值。
而三次样条插值作为一种常用的插值方法，在Matlab中有着广泛的应用。
本文将从简单到复杂，由浅入深地解析Matlab中的三次样条插值方法，以便读者更深入地理解这一技术。
2. 三次样条插值概述三次样条插值是一种利用分段三次多项式对数据点进行插值的方法。
在Matlab中，可以使用spline函数来进行三次样条插值。
该函数需要输入数据点的x和y坐标，然后可以根据需要进行插值操作。
3. 三次样条插值的基本原理在进行三次样条插值时，首先需要对数据点进行分段处理，然后在每个分段上构造出一个三次多项式函数。
这些多项式函数需要满足一定的插值条件，如在数据点处函数值相等、一阶导数相等等。
通过这些条件，可以得到一个关于数据点的插值函数。
4. Matlab中的三次样条插值实现在Matlab中，可以使用spline函数来进行三次样条插值。
通过传入数据点的x和y坐标，可以得到一个关于x的插值函数。
spline函数也支持在已知插值函数上进行插值点的求值，这为用户提供了极大的灵活性。
5. 三次样条插值的适用范围和局限性虽然三次样条插值在许多情况下都能够得到较好的插值效果，但也存在一些局限性。
在数据点分布不均匀或有较大噪音的情况下，三次样条插值可能会出现较大的误差。
在实际应用中，需要根据具体情况选择合适的插值方法。
6. 个人观点和总结通过对Matlab中三次样条插值的深度解析，我深刻地理解了这一插值方法的原理和实现方式。
在实际工程应用中，我会根据数据点的情况选择合适的插值方法，以确保得到准确且可靠的结果。
我也意识到插值方法的局限性，这为我在实际工作中的决策提供了重要的参考。
通过以上深度解析，相信读者已经对Matlab中的三次样条插值有了更加全面、深刻和灵活的理解。
在实际应用中，希望读者能够根据具体情况选择合适的插值方法，以提高工作效率和准确性。
三次样条函数三弯矩算法
摘要求函数在给定区间上的定积分，在微积分学中已给出了许多计算方法，但是，在实际问题计算中，往往仅给出函数在一些离散点的值，它的解析表达式没有明显的给出，或者，虽然给出解析表达式，但却很难求得其原函数。
这时我们可以通过数值方法求出函数积分的近似值。
在用近似值代替真实值时，遇到的问题就是近似值的代数精度是否足够。
当代数精度不足够时，很显然提高插值函数的次数是一种方法，但是考虑到数值计算的稳定性，当次数过高时，会出现龙格现象，用增大n的方法来提高数值积代数精度是不可取的。
正如我们所知道的分段线性插值，逼近程度好，但光滑性差。
分段三次Hermite插值，逼近程度好，光滑性也有所提高，但是需要增加更多的条件，不太实用。
因此，我们将介绍一种结合二者的优点的插值方法——三次样条插值。
本实验将介绍三次样条插值的三弯矩算法。
关键词：龙格现象三弯矩算法代数精度分段三次Hermite插值1、实验目的1) 通过本次实验体会并学习三次样条插值的优点。
2) 通过对三次样条插值进行编程实现，提高自己的编程能力。
3) 用实验报告的形式展现，提高自己在写论文方面的能力。
2、算法流程如果已知函数)(x f y =在节点a =x 0<x 1<⋯<x n =b,y i =f (x i ),i =0,1,2,⋯,n 处的函数值和导数值：n i x f y i i ,,2,1,0),( ==如果)(x S 满足条件：1) S (x )是一个分段的三次多项式且i i y x s =)(；2) S (x )在[a,b]具有二阶连续导数。
则称S (x )是三次样条插值函数。
S (x )的具体形式为：其中S i (x)在[x n−1,x n ]上是三次多项式S i (x )=a i x 3+b i x 2+c i x +d i 由插值条件S (x i )=y i ，i=0,1,2,…,n ，得n+1个条件。
边界条件一：S ′(x 0)=y 0′,S ′(x n )=y n ′边界条件二：S ′′(x 0)=y 0′′,S ′′(x n )=y n ′′边界条件三：假定函数y =f(x)是以b-a 为周期的周期函数，这时要求S(x)也是周期函数，即{S (x 0+0)=S(x n −0)S ′(x 0+0)=S ′(x n −0)S ′′(x 0+0)=S ′′(x n +0)()()()()⎪⎪⎩⎪⎪⎨⎧∈∈∈=-]12,121,01,[,...............][,][,n n n x x x x s x x x x s x x x x s x s针对三种边界条件的求解方法的不同，可以分为三转角算法和三弯矩算法，本实验将介绍和学习三转角算法。
用MATLAB计算等距三次样条插值问题
2 表达式中系数的求解
S 4( π ) 中的任意一个三次样条函数可以表示成
38
n1
四川工业学院学报 2003 年 x ), x ∑ k iB i( ∈ [ a , b] ( 2) 于是求满足条件( 3) 、 ( 4) 的 三次插值样条函数( 2)的 问题转换为求解线性方程组( 7) 的问题 。 只要从( 7)中 解出 k i( i =-1 , 0 , …, n -3) , 即可求得样条函数 。
T
k n -1 = y n 及中间系数满足的等式 k -1 B -1( x 1)+ k 0 B 0( x 1)= y 1 - y 0 + h y′ 0 Bx 1) 2( 3
ki 3 B i3( x i) +k i 2 B i2( xi ) +k i 1 B i1 ( xi )= y i i = 2 , 3 , … , n -2 k n -4 B n -4( xn -1)+k n -3 B n -3 = y n -1 h - y n - y ′ B ( x )= y i 3 n n -2 n -1 ( 6) 利用基函数( 1) , 及已知数据( 3) , 可将( 6) 式写成矩阵 形式 : 7 2 1 4 0 1 1 4 1 1 4 2 1 7 · k -1 k0 k1 ┇ k n -4 k n -3
用matlab计算等距三次样条插值问题matlab等距节点插值三次样条插值matlabmatlab样条插值matlab样条插值函数matlab样条插值求曲率matlabb样条插值拟合matlab中三次样条插值matlabb样条插值双三次样条插值matlab
四川工业学院学报
Journa l of Sichua n University o f Science and Technolog y
计算方法上机作业——求三次样条插值函数的matlab程序
附录 3 求三次样条插值函数的 matlab 程序 for f = 2:n-1; ly = 0; for g = 1:f-1 ly = ly+l(f,g)*yy(g); end yy(f) = D(f)-ly; end M1(n-1) = yy(n-1)/u(n-1,n-1); for rr=1:n-2 r = n-1-rr; uM1 = 0; for s=r+1:n-1 uM1 = uM1+u(r,s)*M1(s); end M1(r) = (yy(r)-uM1)/u(r,r); end M = [M1(n-1,1);M1]; end ss = 0; for t=1:n-1 S(t,1) = (M(t+1)-M(t))/(6*h(t)); S(t,2) = (M(t)*x(t+1)-M(t+1)*x(t))/(2*h(t)); S(t,3) = (M(t+1)*x(t)^2-M(t)*x(t+1)^2)/(2*h(t))+(y(t+1)-y(t))/h(t)+h(t)*(M(t)-M(t+1))/6; S(t,4) = (M(t)*x(t+1)^3-M(1)*x(t)^3)/(6*h(t))+(y(t)*x(t+1)-y(t+1)*x(t))/h(t)+h(t)*(M(t+1)* x(t)-M(t)*x(t+1))/6; for x1 = x(t):(x(t+1)-x(t))/100:x(t+1) ss = ss+1; xx(ss) = x1; SS(ss) = S(t,1)*x1^3+S(t,2)*x1^2+S(t,3)*x1+S(t,4); end end plot(xx,SS,'-k','linewidth',2); hold on plot(x,y,'*k','markersize',10); hold on xlabel('x'); ylabel('S(x)'); grid; fprintf('\n 所求的三次样条插值函数为：\n'); for uu=1:n-1 fprintf('S(x) = %10.5f*x^3+%10.5f*x^2+%10.5f*x+%10.5f, %8.4f<= x <=%8.4f\n',S(uu,1),S(uu,2),S(uu,3),S(uu,4),x(uu),x(uu+1)); end
三次样条插值函数MATLAB编程实现
三次样条插值函数为()()[)()[]1011,,,,n n n S x x x x S x S x x x x-⎧∈⎪=⎨⎪∈⎩ 利用三次埃尔米特插值函数表示三次样条插值函数，即()()()()())111111,,j j j j j j j j j j j S x y x y x m x m x x x x ααββ++++++⎡=+++∈⎣（0,1,,1j n =-）基函数满足()()()()()()21112111121121111212jj j j j j j j j j j j j j j j j j j jj j j j x x x x x x x x xx xx x x x x x xx xx x x x xx x x x x x xααββ++++++++++++⎛⎫⎛⎫--=+ ⎪⎪ ⎪⎪--⎝⎭⎝⎭⎛⎫⎛⎫--=+ ⎪⎪ ⎪⎪--⎝⎭⎝⎭⎛⎫-=-⎪ ⎪-⎝⎭⎛⎫-=-⎪ ⎪-⎝⎭由上式易得()()()()()()()()()()()()()()1331111331112211112211612612246246j j j j j j j j j j j j j j j j j j j j j j jj j j j j x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x xx ααββ+++++++++++++++''=---+''=-+--+''=---+''=---则有()()()()()()()()()()()111111113333111111122221111661212242466j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j S x y x y x m x m x x x x x y x y x x x x x x x x x x x x x m x m x x x x x x x x x ααββ+++++++++++++++++++''''''''''=+++⎡⎤⎡⎤++⎢⎥⎢⎥=-+-+⎢⎥⎢⎥----⎣⎦⎣⎦⎡⎤++⎢⎥+-+-⎢⎥----⎣⎦)1,j j x x x +⎡⎤⎢⎥⎡∈⎣⎢⎥⎣⎦（0,1,,1j n =-）同理有()()()()()()()()()()()()()()()11111113333111111122221111661212242466j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j S x y x y x m x m x x x x x y x y x x x x x x x x x x x x x m x m x x x x x x x x x ααββ------------------''''''''''=+++⎡⎤⎡⎤++⎢⎥⎢⎥=-+-+⎢⎥⎢⎥----⎣⎦⎣⎦⎡⎤⎡++⎢⎥+-+-⎢⎥----⎣⎦⎣)1,j j x x x -⎤⎢⎥⎡∈⎣⎢⎥⎦（1,,j n =）根据样条函数二阶导数连续性，即()()100j j j j S x S x +''''+=-（1,,1j n =-）即()()()()()()()()()()()()()()()()111111332211111111113322111166426624j j jj j j j j j j jj j jj j j j j j jj j j j jj j j j j jjj jj jj jj x x y x x y x x x x m m x x xx xx xx x x y x x y x x x x m m x x xx xx xx ++++++++++--------------+++--------=+++----（1,,1j n =-）化简得()()()()()111111111111233j j j j j j j j j j j j j j j j j j jj j xx m x x m x x m x x x x y y y y x x x x +-+--+-++-+--+-+---=-+---（1,,1j n =-）可得线性方程组()()()()()()()()()()0121201023231213121221111110212110211032213221322122233333n n n n n n n n n n n n m m x x x x x x m x x x x x x m x x x x x x m m m x x x x y y y y x x x x x x x x y y y y x x x x y ------⨯+-+⨯⎛⎫ ⎪ ⎪ ⎪---⎛⎫⎪ ⎪--- ⎪ ⎪⎪ ⎪⎪ ⎪ ⎪--- ⎪⎝⎭ ⎪ ⎪ ⎪⎝⎭---+------+---=()()()121112112113n n n n n n n n n n n n n x x x x y y y x x x x ----------⨯⎛⎫⎪ ⎪ ⎪ ⎪⎪⎪ ⎪-- ⎪-+- ⎪--⎝⎭为了使样条插值问题有惟一解，我们在原有方程基础上增加两个边界条件。
三次样条插值函数的构造与matlab实现
三次样条插值函数的构造与matlab实现
三次样条插值函数是一种常用的数值计算方法，用于在已知一些散点数据的情况下，生成一个平滑的函数来近似这些数据。
其构造方法如下：
1.假设有n个数据点(x1,y1),(x2,y2),...,(xn,yn)，并且这些数据点按照x从小到大的顺序排列。
2.为了使得插值函数在两个相邻的数据点之间是光滑的，需要在每个相邻的数据点(xi,yi)和(xi+1,yi+1)之间插入一个三次多项式Si(x)。
这个三次多项式可以写成以下形式：
Si(x) = ai + bi(x - xi) + ci(x - xi)2 + di(x - xi)3，其中
i = 1, 2, ..., n-1。
3.为了保证插值函数在整个数据区间内光滑，需要添加一些附加条件。
最常用的条件是，插值函数在每个数据点上的一阶导数和二阶导数与相邻插值函数相等。
这些条件可以得到2(n-1)个方程式。
4.解这些方程，得到系数ai, bi, ci, di。
5.最后，插值函数可以表示为：S(x) = Si(x)，其中xi ≤x ≤xi+1。
在Matlab中，可以使用spline函数实现三次样条插值。
LAB03_三次样条插值实验
LAB03_三次样条插值实验实验内容1.用三次样条插值的三弯矩法,编制第一与第二种边界条件的程序.已知数据如下：求)(x f 的三次样条插值函数)(x S 满足： (1)自然边界条件''(0.2)''(1.0)0;S S ==(2)第一种边界条件.55741.1)0.1(',20271.0)2.0('==S S 要求输出用追赶法解出的弯矩向量(0M ,1M ,2M ,3M ,4M )及)1.02.0(i S +(i=0,1,2,3,4,5,6,7,8)的值.并画出)(x S y =的图形.实验原理简述我们先利用一阶导数)(x s ''在节点i x （i=1,2,…n-1）上的连 性以及边界条件，列出确定二阶导数)(i i x s M ''=（i=0,1,2,…n-1）的线性方程组（三弯矩方程），并由此解出i M 然后用i M 表达是)`(x s 即可i i i i i i g M M M =+++-112λμ（三弯矩方程） jj j i h h h +=--11μ，jj j i h h h +=-1λ， ,,,1,0n i =[][][]11111,,6,,6+---+=+-=j j j jj jj j j i x x xf h h x x f x x fg .其中有（1+n ）个未知数n M M M ,...,10，而方程只有（n-1）个，当满足第一种边界条件时，可的另两个方程[]()010010,62f x x f h M M '-=+， []()n n n n n n x x f f h M M ,62111----'=+如果令[]()[]()n n n n n n x x f f h g f x x f h d ,6,1,,6,111010000---'=='-==μλ，将上述方程综合后的一下矩阵形式：222211110nn n μλμλμλ--nn M M M M 110- =n n g g g g 110-可以证明此方程组满足追赶法的条件，我们用追赶法可得M 的)1(+n 值，将其带入公式即得)(x s对第二种边界条件，直接的端点方程n n f M f M ''=''=,00并且令n n n f d f d ''=''===2,2,0000μλ，则又得三弯矩方程同理即可求得解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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;
d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1);
newd=[d0;d];
%%%%%%%%%%%%
function intersanwj(x,y,x0,y0,y1a,y1b)
% 三弯矩样条插值
%第一部分
n=length(x);m=length(y);
if m~=n
error('x or y 输入有误，再来');
end
%重新定义h
h=zeros(n,1);
h(1)=x(1)-x0;
for k=2:n
h(k)=x(k)-x(k-1);
end
%sptep1 调用三弯矩函数
[a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b);
% 三对角方程
M=chase(a,b,c,d);
% 求插值函数
fprintf('三次样条(三弯矩)插值的函数表达式\n');
syms X ;
fprintf('S0--1:\n');
S(1)=collect(((1/6)*M(2)*(X-x0).^3-(1/6)*M(1)*(X-x(1)).^3+(y(1)-
(M(2)*h(1).^2)/6)*(X-x0)-(y0-(M(1)*h(1).^2)/6)*(X-x(1)))/h(1));
for k=2:n
fprintf('S%d--%d:\n',k-1,k);
S(k)=collect(((1/6)*M(k+1)*(X-x(k-1)).^3-(1/6)*M(k)*(X-x(k)).^3+(y(k)-(M(k+1)*h(k).^2)/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k)); end
S=S.';
disp(S);
fprintf('以上为样条函数(三弯矩)解析式，显示为手写如下：\n');
pretty(S);
%第二部分
%是否继续运行程序
myloop=input('继续运行程序输入“1”，否则输入“0”\n');
if myloop
while myloop
xi=input('输入需要计算的点的值，并按回车键\n');
if xi>x0|xi<x(n)
fprintf('现在开始计算输入点的插值函数值……\n');
else
fprintf('输入数值不在插值范围内，请重新输入\n');
xi=input('输入需要计算的点的值，并按回车键……\n');
end
% 确定输入的数值应该使用哪个解析式
newx=[x0;x];
[r,suoy]=min(abs(newx-xi));
fprintf('输入点的插值函数值为:\n\n');fprintf('\t');
if xi<=newx(suoy)
f=subs(S(suoy-1),X,xi);
else
f=subs(S(suoy),X,xi);
end
disp(f);
myloop=input('继续计算输入“1”，终止计算输入“0”\n');
end
else
return;
end
%%%%%%%%%%%%
function [x]=chase(a,b,c,d)
%追赶法解性方程组 a是下三角b是对角线c是上三角 d是常数项%输入的a b c d 均为列向量
n=length(b);
u=zeros(n,1);
v=zeros(n,1);
x=zeros(n,1);
%追
v(1)=c(1)/b(1);u(1)=d(1)/b(1);
for i=2:n-1
v(i)=c(i)/(b(i)-v(i-1)*a(i-1));
u(i)=(d(i)-u(i-1)*a(i-1))/(b(i)-v(i-1)*a(i-1));
end
u(n)=(d(n)-u(n-1)*a(n-1))/(b(n)-v(n-1)*a(n-1)); %赶
x(n)=u(n);
for i=n-1:-1:1
x(i)=u(i)-v(i)*x(i+1);
end。