*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2)+16.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).^2.*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^4.*(2.*cos(2.*x)./(3.+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)))+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^3.*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^4.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3.+4.*x.^2).^2.*x).^4-8.*tan(cos((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2))).*(1.+tan(cos((3.^(1./2)+sin(2.*x))./(3.+4*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3.+4.*x.^2)).^2.*(2.*cos(2.*x)./(3.+4.*x.^2)-8*(3.^(1./2)+sin(2.*x))./(3.+4*x.^2).^2.*x).^4+6.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4-3.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2).^2-4.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)-(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(16.*sin(2.*x)./(3+4.*x.^2)+256.*cos(2.*x)./(3+4.*x.^2).^2.*x-3072.*sin(2.*x)./(3+4.*x.^2).^3.*x.^2+192.*sin(2.*x)./(3+4.*x.^2).^2-24576.*cos(2.*x)./(3+4.*x.^2).^4.*x.^3+3072.*cos(2.*x)./(3+4.*x.^2).^3.*x+98304.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^5.*x.^4-18432.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^2+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3)-12.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).^2*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^3.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)-24.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2.*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))+36.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^2.*cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2)+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*cos((3.^(1./2)+sin(2.*x))./(3
86．
+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).^4+6.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(-4.*sin(2.*x)./(3+4.*x.^2)-32.*cos(2.*x)./(3+4.*x.^2).^2.*x+128.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x.^2-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2).^2+8.*tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).*(1+tan(cos((3.^(1./2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3.^(1./2)+sin(2.*x))./(3+4.*x.^2)).^2.*(2.*cos(2.*x)./(3+4.*x.^2)-8.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^2.*x).*(-8.*cos(2.*x)./(3+4.*x.^2)+96.*sin(2.*x)./(3+4.*x.^2).^2.*x+768.*cos(2.*x)./(3+4.*x.^2).^3.*x.^2-48.*cos(2.*x)./(3+4.*x.^2).^2-3072.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^4.*x.^3+384.*(3.^(1./2)+sin(2.*x))./(3+4.*x.^2).^3.*x)
myxx=max(yxxxx), R=(h.^4).* abs(myxx./384)
将其保存为名为myxx.m的M文件，然后在MATLAB工作窗口输入该文件名
>> myxx 运行后屏幕显示myxx =maxf???x??(4)(x)和Hn,3(x)在区间[??,?]上的误差公式
h4R?maxf384???x??(4)(x)如下
myxx = R =
73.94706841647552 1734520780029061/9007199254740992*h^4
最后在MATLAB工作窗口输入
>>h=2*pi/9; R =1734520780029061/9007199254740992*h^4
运行后屏幕显示Hn,3(x)在区间[??,?]上的误差限
R =
0.04574453029948
2.5 三次样条（Spline）及其MATLAB程序
2.5.4 用一阶导数计算的几种样条函数和MATLAB程序
计算压紧样条的MATLAB主程序
function
[m,H,lambda,mu,D,A,dY,sk]=ClampedSp (X,Y,dyo,dyn) m=length(X);A=zeros(m,m);
n=m-1;H=zeros(1,n);lambda=zeros(1,n);
mu=zeros(1,n);lambda(1)=0;A(1,1)=2;A(1,2)=lambda(1); D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1-lambda(1);D(1)=2*dyo
;
for k=1:n
hk=X(k+1)-X(k);H(k+1)=hk; end
H=H(2:n+1); for k=1:n-1
lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak; muk=1-lambda(k+1);mu(k)=muk;
dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)
-Y(k+1))./H(k+1)));
D(k+1)=dk; end
D(m)=2*dyn;mu(n)=0; n;H;lambda;mu;D; for i=1:m-1
A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); end
dY=A\\D';
87．
syms x m=length(X);S=zeros(m-1,1); for k=2:m
sk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+
Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H(k-1)^2)
end
例2.5.2 用计算压紧样条的MATLAB主程序ClampedSp.m求例6.7.1的问题. 解 保存ClampedSp.m为M文件，输入程序
>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7, [m,H,lambda,mu,D,A,dY,sk]=ClampedSp(X,Y,dyo,dyn)
运行后输出压紧样条函数sk，X的维数m，由hi?1(i?1,2,?,n)、?i、?i和di （6.90）式的系数矩阵A，(i?1,2,?,n?1)的坐标组成的向量H、lambda和mu、D，
线性方程组（6.90）的解向量dY如下
sk =
2*(6-2*x)*x^2+2*x*(x-2)^2+9/4*(x-2)*x^2 sk =
2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-4)^
2+5/2*(x-4)*(x-2)^2
sk =
9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+5/2*(x-4)*(x-6)^2+2*(x-6)*(x-4)^2
sk =
27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x-6)
*(x-10)^2+7/16*(x-10)*(x-6)^2
m = 5 H =
2 2 2 4 lambda =
0 0.5000 0.5000 0.3333 mu =
0.5000 0.5000 0.6667 0 D =
16.0000 27.0000 28.5000 25.0000 14.0000 A =
2.0000 0 0 0 0 0.5000 2.0000 0.5000 0 0 0 0.5000 2.0000 0.5000 0 0 0 0.6667 2.0000 0.3333 0 0 0 0 2.0000 dY =
8.0000 9.0000 10.0000 8.0000 7.0000 输入程序
>> x2=3,
s2=2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+9/4*(x-2)*(x-4)^2+5/2*(x-4)*(x-2)^2
x4=8,
s4=27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+1/2*(x
-6)*(x-10)^2+7/16*(x-10)*(x-6)^2
运行后输出f(3)?s2(3)和f(8)?s4(8)的值如下
88．
x2 = s2 = x4 = s4 = 3 25.7500 8 68.5000 例2.5.2和例2.5.1计算的结果完全相同. 计算自然样条的MATLAB主程序
function [m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y) m=length(X);A=zeros(m,m);
n=m-1;H=zeros(1,n);lambda=zeros(1,n);
mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1); D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1)); for k=1:n
hk=X(k+1)-X(k);H(k+1)=hk; end
H=H(2:n+1); for k=1:n-1
lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak; muk=1-lambda(k+1);mu(k)=muk;
dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)
-Y(k+1))./H(k+1)));
D(k+1)=dk; end
D(m)= 3*(Y(m)-Y(m-1))/H(m-1);mu(n)=1; n;H;lambda;mu;D; for i=1:m-1
A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); end
dY=A\\D'; syms x
m=length(X);S=zeros(m-1,1); for k=2:m
sk=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+
Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H(k-1)^2)
end
例2.5.3 用计算自然样条的MATLAB主程序NaturalSp.m求例6.7.1的问题. 解 保存ClampedSp.m为M文件，输入程序
>> X=[0:2:6,10],Y=[0,16,36,54,82],dyo=8,dyn=7, [m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)
运行后输出自然样条函数sk， X的维数m，由hi?1(i?1,2,?,n)、?i、?i和di （6.88）式的系数矩阵A，(i?1,2,?,n?1)的坐标组成的向量H、lambda和mu、D，
线性方程组（6.88）的解向量dY如下
sk =
2*(6-2*x)*x^2+915/172*x*(x-2)^2+117/86*(x-2)*x^2 sk =
2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+117/86*(x-2)*(x-4)^2+471/172*(x-4)*(x-2)^2
sk =
9/2*(-6+2*x)*(x-6)^2+27/4*(14-2*x)*(x-4)^2+471/172*(x-4)
*(x-6)^2+333/172*(x-6)*(x-4)^2
sk =
27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+333/688*(
x-6)*(x-10)^2+285/688*(x-10)*(x-6)^2
m = 5
H = 2 2 2 4
lambda = 1 1/2 1/2 1/3 mu = 1/2 1/2 2/3 1
D = 48 27 57/2 25 21 A =
89．
2 1 0 0 0 1/2 2 1/2 0 0 0 1/2 2 1/2 0 0 0 2/3 2 1/3 0 0 0 1 2 dY =
915/43 234/43 471/43 333/43 285/43 输入程序
>> x2=3,
s2=2*(-2+2*x)*(x-4)^2+9/2*(10-2*x)*(x-2)^2+117/86*(x-2)*
(x-4)^2+471/172*(x-4)*(x-2)^2
x4=8,
s4=s4=27/32*(-8+2*x)*(x-10)^2+41/32*(24-2*x)*(x-6)^2+333
/688*(x-6)*(x-10)^2+285/688*(x-10)*(x-6)^2 运行后输出f(3)?s2(3)和f(8)?s4(8)的值如下
x2 = s2 = x 4= s4 = 3 24.6221 8 68.5581
例2.5.3和例2.5.2用的方法不同，所以计算的结果不同，但是两种方法计算的结果相近.
2.5.6 用MATLAB计算三次样条
例2.5.6 给出节点的数据如下： -2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31 分别求在下列条件下在插值点x1??0.02，x2?2.56处的压紧三次样条插值，并显示该样条函数的有关信息：
?(?1)?5，Sn?(3.50)?29.16； （1）端点约束条件为Snx -1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50 f(x) ?(?1)?0，Sn?(3.50)?0. （2）端点约束条件为Sn解 （1）输入MATLAB程序
>> X=[-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82
3.50];
Y=[-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31];
XI=[-0.02 2.56]; YI= spline (X, [5,Y,29.16],XI), PP = spline (X, [5,Y,29.16])
x2?2.56处的插值和该样条函数的有关信运行后屏幕显示压紧样条分别在x1??0.02，
息如下
YI =
-3.1058 4.7834 PP = form: 'pp'
breaks: [-1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000] coefs: [8x4 double] pieces: 8
90．