计算方法大作业（第二次）(2)

A=zeros(n);B=zeros(n,1); for i=1:n-1 a(i)=y(i);

dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i); end

for i=2:n-1

A(i,i-1)=dx(i-1);

A(i,i)=2*(dx(i-1)+dx(i)); A(i,i+1)=dx(i);

B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1)); end

%自然样条端点条件（端点二阶导数为0） if flag==0 A(1,1)=1; A(n,n)=1; end

%端点一阶导数条件 if flag==1

A(1,1)=2*dx(1);A(1,2)=dx(1);

A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1); B(1,1)=3*(dy(1)/dx(1)-vl); B(n,1)=3*(vr-dy(n-1)/dx(n-1)); end

%端点二阶导数条件 if flag==2

A(1,1)=2;A(n,n)=2; B(1,1)=vl;B(n,1)=vr; end c=A\\B; for(i=1:n-1)

d(i)=(c(i+1)-c(i))/(3*dx(i));

b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3; end

[mm,nn]=size(xx); yy=zeros(mm,nn); for i=1:mm*nn for ii=1:n-1

if xx(i)>=x(ii)&xx(i)

elseif xx(i)==x(n) j=n-1; end

end

yy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3; end end