李庆扬数值分析MATLAB代码

function yy = spline( x, z, flag,xx ) % 三次样条插值，x为列数据,z为端点导数值 %flag=1,已知两端导数值;flag=2,已知两端二阶导数；flag=3,周期函数 [n, m] = size(x); a = x(:,1); b = x(:,2); for j=1:(n-1) h(j) = a(j+1)-a(j); end for j=1:n-2 u(j) = h(j)/(h(j)+h(j+1)); end for j = 2:(n-1) v(j) = h(j)/(h(j-1)+h(j)); f1(j)= (b(j+1)-b(j))/(a(j+1)-a(j)); f2(j)= (b(j)-b(j-1))/(a(j)-a(j-1)); d(j) = 6*(f1(j)-f2(j))/(h(j-1)+h(j)); end switch flag case 1 A(n,n) = 0; v(1) = 1; u(n-1) = 1; d(1) = 6*((b(2)-b(1))/(a(2)-a(1))-z(1))/h(1); d(n) = 6*(z(2)-(b(n)-b(n-1))/(a(n)-a(n-1)))/h(n-1); A(1,:) = [2,v(1),zeros(1,n-2)]; for i=2:(n-1) A(i,:)=[zeros(1,i-2),u(i-1),2,v(i),zeros(1,n-i-1)]; end A(n,:)=[zeros(1,n-2),u(n-1),2]; case 2 A(n,n) = 0; v(1) = 0; u(n-1)= 0; d(1) = 2*z(1); d(n) = 2*z(2); A(1,:) = [2,v(1),zeros(1,n-2)]; for i=2:(n-1) A(i,:)=[zeros(1,i-2),u(i-1),2,v(i),zeros(1,n-i-1)]; end A(n,:)=[zeros(1,n-2),u(n-1),2]; case 3 v(n)= h(1)/(h(1)+h(n-1)); u(n-1)= 1-v(n); f1 = (b(2)-b(1))/(a(2)-a(1)); f2 = (b(n)-b(n-1))/(a(n)-a(n-1)); d(n)= 6*(f1-f2)/(h(1)+h(n-1)); A(1,:) = [2,v(2),zeros(1,n-4),u(1)]; for i=2:(n-2) A(i,:)=[zeros(1,i-2),u(i),2,v(i+1),zeros(1,n-i-2)]; end A(n-1,:)=[v(n),zeros(1,n-4),u(n-1),2]; d = d(2:n); end D = d'; M = linsolve(A,D); if flag==3 M(1)=M(n-1); M(2:n)=M(1:n-1); end % syms t S % S = []; for i =1:n-1 s(i,1) = M(i)/6/h(i); s(i,2) = (b(i)-M(i)*h(i)^2/6)/h(i); s(i,3) = M(i+1)/6/h(i); s(i,4) = (b(i+1)-M(i+1)*h(i)^2/6)/h(i); % clear S % S(i) = s(i,1)*(a(i+1)-t)^3 + s(i,2)*(a(i+1)-t) + s(i,3)*(t-a(i))^3 + s(i,4)*(t-a(i)); end figure for i=1:n-1 if (a(i)<=xx && xx<=a(i+1)) yy = s(i,1)*(a(i+1)-xx)^3 + s(i,2)*(a(i+1)-xx) +s(i,3)*(xx-a(i))^3 + s(i,4)*(xx-a(i)); end xi = a(i):0.01:a(i+1); for j=1:length(xi) yi(j) = s(i,1)*(a(i+1)-xi(j))^3 + s(i,2)*(a(i+1)-xi(j)) + s(i,3)*(xi(j)-a(i))^3 + s(i,4)*(xi(j)-a(i)); ci(i) = min(yi); di(i) = max(yi); end c = min(ci); d = max(di); plot(xi,yi) hold on axis([a(1)-(a(n)-a(1))/13 a(n)+(a(n)-a(1))/13 c-(d-c)/13 d+(d-c)/13]); end scatter(a,b,'r*') hold on scatter(xx,yy,'b') end