cubic_line1.m

数值分析

function yi=cubic_line1(x,y,ydot,xi)
n=length(x);ny=length(y);m=length(xi);
if n~= ny
    error('The lengths of Xand Y must be equal ');
    return;
end
if isempty(ydot)==1
    ydot=[(y(2)-y(1))/(x(2)-x(1)) (y(n)-y(n-1))/(x(n)-x(n-1))];
end
h=zeros(1,n);lambda=ones(1,n);mu=ones(1,n);
M=zeros(n,1);d=zeros(n,1);
for k=2:n
    h(k)=x(k)-x(k-1);
    if abs(h(k))<eps
        error('the DATA is error!');
        return;
    end
end
for k=2:n-1
    lambda(k)=h(k+1)/(h(k)+h(k+1));
    mu(k)=1-lambda(k);
    d(k)=6/(h(k)+h(k+1))*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/h(k));
end
d(1)=6/h(2)*((y(2)-y(1))/h(2)-ydot(1));
d(n)=6/h(n)*(ydot(2)-(y(n)-y(n-1))/h(n));
A=diag(2*ones(1,n));
for i=1:n-1
    A(i,i+1)=lambda(i);A(i+1,i)=mu(i+1);
end
M=A\d;
for k=2:n
    for j=1:m
        if x(k-1)<=xi(j) & xi(j)<=x(k)
        yi(j) = M(k-1)/6/h(k)*(x(k)-xi(j))^3+M(k)/6/h(k)*(xi(j)-x(k-1))^3+1/h(k)*(y(k)-M(k)*h(k)^2/6)*(xi(j)-x(k-1))+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-xi(j));
        return;
    end
end
end