第3节 函数近似方法(附录).docVIP

第3章 附录 3.1.5 n+1点n次插值 function y = lagrange_2(x0,y0,x) n = length(x0); m = length(x); for i = 1:m z = x(i); s = 0.0; for k = 1:n p = 1.0; for j =1:n if j~=k p = p*(z-x0(j))/(x0(k)-x0(j)); end end s = s+ p*y0(k); end y(i) = s; End ============================================== 3.1.6 三次样条插值 例题3.1.7计算程序 % demo_spline.m x = -1:0.1:1; y = 1./(1+x.*x); xx = -1:.01:1; yy = spline(x,y,xx); plot(x,y,bo,xx,yy,r,LineWidth,2) title(y=1/(1+x^2): 样条插值… ,FontSize,12); xlabel(x,FontSize,12); ylabel(y,FontSize,12) % demo_splines.m x = [-1.:0.1:1.]; n = length(x); y = 1./(1.+x.*x); m = 10*(n-1)+1; dx = (x(n)-x(1))/(m-1.) for i = 2:n-1 a(i) = x(i+1)-x(i); b(i) = 2.*(x(i+1)-x(i)); c(i) = x(i+1)-x(i); f(i) = 6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i+1)-x(i))); end a(1) = 0.; b(1) = 1.;c(1) = 0.; f(1) = 0.; a(n) = 0.; b(n) = 1.;c(n) = 0.; f(n) = 0.; ypp = tri(a,b,c,f); for i = 1:n-1 c1(i) = (ypp(i+1)-ypp(i))/(x(i+1)-x(i)); c2(i) = (x(i+1)*ypp(i)-x(i)*ypp(i+1))/(x(i+1)-x(i)); c3(i) = ((x(i)*x(i)-2.*x(i+1)*x(i+1)-2.*x(i)*x(i+1))*ypp(i)+(2.*x(i)*x(i)-x(i+1)*x(i+1)... + 2.*x(i)*x(i+1))*ypp(i+1)+6.*(y(i+1)-y(i)))/6./(x(i+1)-x(i)); c4(i) = -(x(i)*x(i+1)*((x(i)-2.*x(i+1))*ypp(i)+(2.*x(i)-x(i+1))*ypp(i+1))... + 6.*(x(i)*y(i+1)-x(i+1)*y(i)))/6./(x(i+1)-x(i)); end j =1; for i = 1:n-1 dj = floor((x(i+1)-x(i)+dx/2)/dx); for k =j:(j+dj-1) z(k) = x(1)+(k-1)*dx; s(k) = c1(i)*z(k)*z(k)*z(k)/6.+c2(i)*z(k)*z(k)/2.+c3(i)*z(k)+c4(i); sp(k)= c1(i)*z(k)*z(k)/2.+c2(i)*z(k)+c3(i); spp(k)=c1(i)*z(k)+c2(i); end j = j+dj; end yy = 1./(1+z.*z); plot(z,yy,k,x,y,bo,z,s,r,LineWidth,2) xlabel(x,FontSize,12); ylabel(y,FontSize,12) legend(原曲线,插值点,插值曲线); title(y=1/(1+x^2):样条插值,FontSize,12); ============================================== 例题3.2.1计算程序 !!!!linear_fit.for!!! program linear_fit ! linear_fit implicit none integer n,i,m parameter (n=7,m=20)