已知一组数据点，编写一程序求解三次样条插值函数满足
并针对下面一组具体实验数据
求解，其中边界条件为.
解：Matlab计算程序为：
clear
clc
x=[0.25 0.3 0.39 0.45 0.53]
y=[0.5000 0.5477 0.6245 0.6708 0.7280]
n=length(x);
for i=1:n-1
h(i)=x(i+1)-x(i);
end
for i=1:n-2
k(i)=h(i+1)/(h(i)+h(i+1));
u(i)=h(i)/(h(i)+h(i+1));
end
for i=1:n-2
gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i));
end
g0=3*(y(2)-y(1))/h(1);
g00=3*(y(n)-y(n-1))/h(n-1);
g=[g0 gl g00];
g=transpose(g)
k1=[k 1];
u1=[1 u];
Q=2*eye(5)+diag(u1,1)+diag(k1,-1)
m=transpose(Q\g)
syms X;
for i=1:n-1
p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i);
p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1);
p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i);
p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1);
p(i)=p1(i)+p2(i)+p3(i)+p4(i);
p(i)=simple(p(i));
end
s1=p(1)
s2=p(2)
s3=p(3)
s4=p(4)
for k=1:4
for z=x(k):0.001:x(k+1)
q=eval(subs(p(k),'X','z'));
plot(z,q,'b')
hold on
end
end
grid on
legend('èy′??ùì??ú??')
title('2??μ')
xlabel('x')
ylabel('p')
计算结果为：
s1 =
-(705394867539368680*X^3 - 529046150654530286*X^2 + 23087199381998953*X - 40023205577025431)/112589990684262400
s2 =
(257361898089296225*X^3)/136796838681378816 - (160081743506404901*X^2)/60798594969501696 + (404209705972252727*X)/202661983231672320 + 429142243010323951/3166593487994880000
s3 =
- (3495912536773825*X^3)/7599824371187712 + (1437374409830143*X^2)/13510798882111488 + (10427488839800859*X)/11258999068426240 + 12358231431982259943/45035996273704960000
s4 =
(38626753769033575*X^3)/18014398509481984 - (245666153971053021*X^2)/72057594037927936 +
(3614707928905781673*X)/1441151880758558720 + 26732501105874704913/720575940379279360000
通过绘制曲线，发现为3次样条曲线，且数据拟合较好。
计算方法实验—三次样条插值 机电学院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] 画图验证：
（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三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（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;
MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤
（第一边界条件）源代码：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)
非常类似前面的三弯矩法，这里的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) % 三转角样条插值
(第一边界条件源代码: 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 (自然边界条件源代码: 仅仅需要对上面部分标注的位置做如下修改 :
(){} 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 ，算法组织如下：
Ｍ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 ，算法组织如下：
告： 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; %检错
x0=[0 0.9211 1.8431 2.9497 3.8714 4.9781 5.9 7.0064 7.9286 8.9678 10.9542 12.0328 12.9544 13.8758 14.9822 15.9039 16.8261 17.9317 19.0375 19.9594 20.8392 22.9581 23.88 24.9869 25.9083]; >> >> y0=[14405 11180 10063 11012 8797 9992 8124 10160 8488 11018 19469 20196 18941 15903 18055 15646 13741 14962 16653 14496 14648 15225 15264 13708 9633]; >> x=0:0.1:25.9; >> y1=interp1(x0,y0,x,'spline'); >> pp1=csape(x0,y0); %样条插值工具箱函数 y2=ppval(pp1,x); %计算x对应的y值 pp2=csape(x0,y0,'second'); y3=ppval(pp2,x); xydata=[x',y1',y2',y3'] subplot(1,2,1) plot(x0,y0,'+',x,y1) title('Spline1') subplot(1,2,2) plot(x0,y0,'+',x,y2) title('Spline2') dx=diff(x); dy=diff(y2); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y2(13<=x&x<=15); ymin=min(ytemp); xmin=x(y2==ymin); xymin_1315=[xmin,ymin]
三次样条插值多项式 ——计算物理实验作业四 陈万物理学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);
4多项式插值与函数最佳逼近 37（上机题）3次样条插值函数： （1）编制求第一型3次样条插值函数的通用程序；（2）已知汽车门曲线型值点的数据如下： 端点条件为8.0' 0=y ，2.0' 10=y ，用所编程序求车门的3次样条插值函数S （x ），并打印出9,,1,0),5.0(?=+i i S 。用matlab 编写 通用程序为： function [Sx ]=Threch(X,Y,dy0,dyn ) %X 为输入变量x 的数值%Y 为函数值y 的数值%dy0为左端一阶导数值%dyn 为右端一阶导数值%Sx 为输出的函数表达式 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);end for i=2:n %求函数的二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1));d(i)=6*f2(i);end d(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-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2;n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; syms x; for i=1:n Sx(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)); end for i=1:n disp('S(x)='); fprintf('%s(%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end S=zeros(1,n); for i=1:n x=X(i)+0.5; S(i)=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; end disp('S(i+0.5)'); disp('i X(i+0.5)S(i+0.5)'); for i=1:n fprintf('%d%.4f%.4f\n',i,X(i)+0.5,S(i)); end End 在运行窗口输入： >>X=[012345678910];Y=[2.513.304.044.705.225.545.785.405.575.705.80]; Threch(X,Y,0.8,0.2)
实验四用MATLAB实现拉格朗日插值、分段线性插值 一、实验目的： 1）学会使用MATLAB软件； 2）会使用MATLAB软件进行拉格朗日插值算法和分段线性差值算法； 二、实验内容： 1用MATLAB实现y = 1./(x.^2+1);（-1<=x<=1）的拉格朗日插值、分段线性 2.选择以下函数，在n个节点上分别用分段线性和三次样条插值的方法，计算m个插值点的函数值，通过数值和图形的输出，将插值结果与精确值进行比较，适当增加n，再作比较，由此作初步分析： (1).y=sinx;( 0≤x≤2π) (2).y=(1-x^2)（-1≤x≤1） 三、实验方法与步骤： 问题一用拉格朗日插值法 1）定义函数：y = 1./(x.^2+1);将其保存在f.m 文件中，程序如下： function y = f1(x) y = 1./(x.^2+1); 2）定义拉格朗日插值函数：将其保存在lagrange.m 文件中，具体实现程序编程如下：function y = lagrange(x0,y0,x) m = length(x); /区间长度/ n = length(x0); for i = 1:n l(i) = 1; end for i = 1:m for j = 1:n for k = 1:n if j == k continue; end l(j) = ( x(i) -x0(k))/( x0(j) - x0(k) )*l(j); end end end y = 0; for i = 1:n y = y0(i) * l(i) + y; end 3）建立测试程序，保存在text.m文件中，实现画图：
已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 求解，其中边界条件为. 解：Matlab计算程序为： clear clc x=[0.25 0.3 0.39 0.45 0.53] y=[0.5000 0.5477 0.6245 0.6708 0.7280] n=length(x); for i=1:n-1 h(i)=x(i+1)-x(i); end for i=1:n-2 k(i)=h(i+1)/(h(i)+h(i+1)); u(i)=h(i)/(h(i)+h(i+1)); end for i=1:n-2 gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end g0=3*(y(2)-y(1))/h(1); g00=3*(y(n)-y(n-1))/h(n-1); g=[g0 gl g00]; g=transpose(g) k1=[k 1]; u1=[1 u]; Q=2*eye(5)+diag(u1,1)+diag(k1,-1) m=transpose(Q\g) syms X; for i=1:n-1 p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i); p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1); p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i); p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1);
并针对下面一组具体实验数据 求解，其中边界条件为. 解：Matlab计算程序为： clear clc x=[0.25 0.3 0.39 0.45 0.53] y=[0.5000 0.5477 0.6245 0.6708 0.7280] n=length(x); for i=1:n-1 h(i)=x(i+1)-x(i); end for i=1:n-2 k(i)=h(i+1)/(h(i)+h(i+1)); u(i)=h(i)/(h(i)+h(i+1)); end for i=1:n-2 gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end g0=3*(y(2)-y(1))/h(1); g00=3*(y(n)-y(n-1))/h(n-1); g=[g0 gl g00]; g=transpose(g) k1=[k 1]; u1=[1 u]; Q=2*eye(5)+diag(u1,1)+diag(k1,-1) m=transpose(Q\g) syms X; for i=1:n-1 p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i); p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1); p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i); p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1); p(i)=p1(i)+p2(i)+p3(i)+p4(i); p(i)=simple(p(i)); end s1=p(1) s2=p(2) s3=p(3) s4=p(4) for k=1:4
23、汽车门曲线三次样条插值曲线相关程序以及结果 原始数据点： x = 0:10; %取自变量为1,2,3， (10) y = [2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80]; %输入因变量y的值 xx = linspace(min(x),max(x),200); %在x的上下界之间取200个插值节点 pp = csape(x,y,'comlete',[0.8,0.2]); %分段三次样条插值，边界条件为左右端点的一阶导数为0.8和0.2 yy = ppval(pp,xx);%计算200个插值节点对应的y值 plot(x,y,'ko',xx,yy,'k') %画出给定的11个点以及插值函数的图像
24、飞鸟外形上部自然边界条件的三次样条插值曲线相关程序以及结果 原始数据如下： x =[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3]; y = [1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];
xx = linspace(min(x),max(x),200); pp = csape(x,y,'second'); %分段三次样条插值，边界条件为左右端点的二阶导数为0，也称为自然边界条件 yy = ppval(pp,xx); plot(x,y,'ko',xx,yy,'k')
clear; f = input('请输入函数表达式:f(x) = ', 's');%注's'，表明允许用户输入一个字符串 a = input('请输入区间左端值a：'); b = input('请输入区间右端值b：'); n = input('请输入区间等分值n：'); for i=1:n+1 x(i) = a + (b-a)/n*(i-1); y(i) = eval(subs(f,'x(i)','x')); end n=n+1; lamda(1)=1;%构造向量 miu(n)=1; %构造向量 h=diff(x); %若X为向量，Y = diff（X）= [X(2)-X(1),X(3)-X(2),...,X(n)-X(n-1)] 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 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)); 它的返回值是符号多项式的系数，依次输出由高阶到0阶的系数，如matlab提供的示例：syms x;sym2poly(x^3 - 2*x - 5)结果为：1 0 -2 -5 end %生成三次样条插值函数图象 for j=1:n-1 x1=x(j):0.01:x(j+1); y1=polyval(S(j,:),x1);多项式的估值运算，使用方法y = polyval(p,x)，返回n 次多项式p在x处的值。输入变量p是一个长度为n+1的向量，其元素为按降幂排列的多项式系数。y=p1*x^n+p2*x^(n-1)+...+pn*x+p(n+1) plot(x1,y1,'b',x,y,'r'); title('spline 三次样条插值函数图象');