初学数值分析——关于三次样条逼近
%三次样条逼近%由x计算h;
x=0:10;
y=;
h=zeros(1,10);
for i=1:10,
h(i)=(x(i+1)-x(i));
end
%计算a,b;
a=zeros(1,11)
b=zeros(1,11)
for i=2:10,
a(i)=(h(i-1)/(h(i-1)+h(i))),
b(i)=3*((1-a(i))*(y(i)-y(i-1))/h(i-1)+a(i)*(y(i+1)-y(i))/h(i))
end
%A
for i=1:9,
A(i,i)=1-a(i)
A(i,i+1)=2
A(i,i+2)=a(i)
end
%m
M=b(2:10)'\A
%s
m=M
t=sym('t')
figure(1);
AXIS()
for i=1:10,
s=(1+2*(t-x(i))/(x(i+1)-x(i)))*((t-x(i+1))/(x(i)-x(i+1)))^2*y(i)+(1-2*(t-x(i+1))/(x(i)-x(i+1)))*((t-x(i))/(x(i+1)-x(i)))^2*y(i+1)+(t-x(i))*((t-x(i+1))/(x(i)-x(i+1)))^2+m(i)+(t-x(i+1))*((t-x(i))/(x(i+1)-x(i)))^2*m(i+1);
ezplot(s,)
hold on;
AXIS()
end
grid on;
这是我按照徐翠薇的《计算方法引论》编得一个matlab程序，可是我在matlab中查spindle函数的帮助，却什莫也看不懂，能帮我解释一下吗？
function output = spline(x,y,xx)
%SPLINE Cubic spline data interpolation.
% YY = SPLINE(X,Y,XX) uses cubic spline interpolation to find YY, the values
% of the underlying function Y at the points in the vector XX.The vector X
% specifies the points at which the data Y is given.If Y is a matrix, then
% the data is taken to be vector-valued and interpolation is performed for
% each column of Y and YY will be length(XX)-by-size(Y,2).
%
% PP = SPLINE(X,Y) returns the piecewise polynomial form of the cubic spline
% interpolant for later use with PPVAL and the spline utility UNMKPP.
%
% Ordinarily, the not-a-knot end conditions are used. However, if Y contains
% two more values than X has entries, then the first and last value in Y are
% used as the endslopes for the cubic spline.Namely:
% f(X) = Y(:,2:end-1), df(min(X)) = Y(:,1), df(max(X)) = Y(:,end)
%
% Example:
% This generates a sine curve, then samples the spline over a finer mesh:
% x = 0:10;y = sin(x);
% xx = 0:.25:10;
% yy = spline(x,y,xx);
% plot(x,y,'o',xx,yy)
%
% Example:
% This illustrates the use of clamped or complete spline interpolation where
% end slopes are prescribed. Zero slopes at the ends of an interpolant to the
% values of a certain distribution are enforced:
% x = -4:4; y = ;
% cs = spline(x,);
% xx = linspace(-4,4,101);
% plot(x,y,'o',xx,ppval(cs,xx),'-');
%
% See also INTERP1, PPVAL, SPLINES (The Spline Toolbox).
% Carl de Boor 7-2-86
% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 5.18 $$Date: 2002/06/06 13:39:51 $
% Generate the cubic spline interpolant in ppform, depending on the
% number of data points (and usually using the not-a-knot end condition).
output=[];
n=length(x);
if n<2, error('There should be at least two data points.'), end
if any(diff(x)<0), =sort(x); else, ind=1:n; end
x=x(:); dx = diff(x);
if all(dx)==0, error('The data abscissae should be distinct.'), end
= size(y); % if Y happens to be a column matrix, change it to
% the expected row matrix.
if yn==1, yn=yd; y=reshape(y,1,yn); yd=1; end
if yn==n
notaknot = 1;
elseif yn==n+2
notaknot = 0; endslopes = y(:,).'; y(:,)=[];
else
error('Abscissa and ordinate vector should be of the same length.')
end
yi=y(:,ind).'; dd = ones(1,yd);
dx = diff(x); divdif = diff(yi)./dx(:,dd);
if n==2
if notaknot, % the interpolant is a straight line
pp=mkpp(x.',,yd);
else % the interpolant is the cubic Hermite polynomial
divdif2 = diff()./dx(,dd);
pp = mkpp(x,...
[(diff(divdif2)./dx(1,dd)).' (*divdif2).' ...
endslopes(1,:).' yi(1,:).'],yd);
end
elseif n==3¬aknot, % the interpolant is a parabola
yi(2:3,:)=divdif;
yi(3,:)=diff(divdif)/(x(3)-x(1));
yi(2,:)=yi(2,:)-yi(3,:)*dx(1);
pp = mkpp(,yi(,:).',yd);
else % set up the sparse, tridiagonal, linear system for the slopes atX .
b=zeros(n,yd);
b(2:n-1,:)=3*(dx(2:n-1,dd).*divdif(1:n-2,:)+dx(1:n-2,dd).*divdif(2:n-1,:));
if notaknot
x31=x(3)-x(1);xn=x(n)-x(n-2);
b(1,:)=((dx(1)+2*x31)*dx(2)*divdif(1,:)+dx(1)^2*divdif(2,:))/x31;
b(n,:)=...
(dx(n-1)^2*divdif(n-2,:)+(2*xn+dx(n-1))*dx(n-2)*divdif(n-1,:))/xn;
else
x31 = 0; xn = 0; b(,:) = dx(,dd).*endslopes;
end
c = spdiags([ ...
;dx(n-2)] ...
],[-1 0 1],n,n);
% sparse linear equation solution for the slopes
mmdflag = spparms('autommd');
spparms('autommd',0);
s=c\b;
spparms('autommd',mmdflag);
% convert to pp form
c4=(s(1:n-1,:)+s(2:n,:)-2*divdif(1:n-1,:))./dx(:,dd);
c3=(divdif(1:n-1,:)-s(1:n-1,:))./dx(:,dd) - c4;
pp=mkpp(x.', ...
reshape([(c4./dx(:,dd)).' c3.' s(1:n-1,:).' yi(1:n-1,:).'], ...
(n-1)*yd,4),yd);
end
if nargin==2
output=pp;
else
output=ppval(pp,xx);
end 原帖由 vib 于 2007-4-4 15:48 发表
%三次样条逼近
%由x计算h;
x=0:10;
y=;
h=zeros(1,10);
for i=1:10,
h(i)=(x(i+1)-x(i));
end
%计算a,b;
a=zeros(1,11)
b=zeros(1,11)
for i=2:10,
a(i)=(h(i-1)/(h(i-1)+h(i))),
b(i) ...
要全面了解matlab的实现，建议阅读 DeBoor 的 A practical guide to Splines 一书 MATLAB里面的函数都经过优化,在原来基础上可能会有不少改变,所以它并不一定按书上写的方法一模一样.
页:
[1]