做面包有哪些网站知乎江北seo页面优化公司

function [y0,C,x]=Spline_interp_1(X,Y,s0,sN,x0)

%X,Y是已知插值点坐标

%s0,sN是两端点的一次导数值

%x0是插值点

%y0是三次样条函数在x0处的值

%C是分段三次样条函数的系数,按升幂排列,区间[Xj,Xj+1]上的表达式为C(j,1)+C(j,2)*(x-Xj)++C(j,3)*(x-Xj)^2+C(j,4)*(x-Xj)^3.

%x是用追赶法依下标顺利求得的插值节点X,Y处的二阶导数值

N=length(X);

C=zeros(4,N-1); h=zeros(1,N-1);

mu=zeros(1,N-1); lmt=zeros(1,N-1);

d=zeros(1,N); %d表示右端函数值

h=X(1,2:N)-X(1,1:N-1);

mu(1,N-1)=1; lmt(1,1)=1;

mu(1,1:N-2)=h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1));

lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1:N-2)+h(1,2:N-1));

d(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);

d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);

d(1,2:N-1)=6*((Y(1,3:N)-Y(1,2:N-1))./h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2))./h(1,1:N-2))./(h(1,1:N-2)+h(1,2:N-1));

%追赶法解三对角方程组

bit=zeros(1,N-1);

bit(1,1)=lmt(1,1)/2;

for i=2:N-1

bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1));

end

y=zeros(1,N);

y(1,1)=d(1,1)/2;

for i=2:N

y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1))/(2-mu(1,i-1)*bit(1,i-1));

end

x=zeros(1,N);

x(1,N)=y(1,N);

for i=N-1:-1:1

x(1,i)=y(1,i)-bit(1,i)*x(1,i+1);

end

v=zeros(1,N-1);

v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1))./h(1,1:N-1);

C(1,:)=Y(1,1:N-1);

C(2,:)=v-h.*(2*x(1,1:N-1)+x(1,2:N))/6;

C(3,:)=x(1,1:N-1)/2;

C(4,:)=(x(1,2:N)-x(1,1:N-1))./(6*h);

if nargin<5

y0=0;

else

for j=1:N-1

if x0>=X(1,j)& x0<=X(1,j+1)

omg=x0-X(1,j);

y0=((C(4,j)*omg+C(3,j))*omg+C(2,j))*omg+C(1,j);

end

end

end

C=C'

取消

评论