第3章函数近似方法(附录

1、第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.mx = -1:0.1:1; y = 1./(1+x.*x);xx = -1:.0

2、1:1;yy = spline(x,y,xx);plot(x,y,bo,xx,yy,r,LineWidth,2)title(y=1/(1+x2): 样条插值,FontSize,12);xlabel(x,FontSize,12);ylabel(y,FontSize,12)% demo_splines.mx = -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) =

3、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);enda(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

4、(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);endj =1;for i = 1:n-1 dj = floor(x(i+1)-x(i)+dx/2)/d

5、x); 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;endyy = 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

6、(原曲线,插值点,插值曲线); title(y=1/(1+x2):样条插值,FontSize,12);=例题3.2.1计算程序!linear_fit.for! program linear_fit! linear_fitimplicit none integer n,i,mparameter (n=7,m=20)real, dimension (n): x,yreal, dimension (m): xi,yireal a,b,dxopen(5,file=exa3_2_1.dat) open(2,file=exa3_2_1_old.dat)data x/0.5, 1.2, 2.1, 2.9,

7、3.6, 4.5, 5.7/data y/2.81,3.24,3.80,4.30,4.73,5.29,6.03/ write(5,(2f14.6) (x(i),y(i),i=1,n) call fit(n,x,y,a,b)print *,a=,a, b=,bdx=(x(n)-x(1)/(m-1)do i=1,m xi(i)=x(1)+(i-1)*dx yi(i)=a+b*xi(i)end do write(2,(2f14.6) (xi(i),yi(i),i=1,m)endsubroutine fit(n,x,y,a,b)implicit noneinteger n,ireal, dimensi

8、on (n): x,y real xp,yp,x2,xyp,a,b,cxp = 0.0yp = 0.0x2 = 0.0xyp= 0.0do i = 1,n x2 = x2+x(i)*x(i) xp = xp+x(i) yp = yp+y(i) xyp= xyp+x(i)*y(i)end do xp = xp/nyp = yp/nx2 = x2/nxyp= xyp/nc = x2-xp*xp a = (yp*x2-xp*xyp)/cb = (xyp-xp*yp)/c return end/*linear fit.cpp*/#include main() int n =7,i; float x=0

9、.5, 1.2, 2.1, 2.9, 3.6, 4.5, 5.7; float y=2.81,3.24,3.80,4.30,4.73,5.29,6.03; float sx=0.,sy=0.,sxy=0.,sx2=0,deno,a,b; for(i=0;i=6;i+) sx=sx+xi; sy=sy+yi; sxy=sxy+xi*yi; sx2=sx2+xi*xi; deno=n*sx2-sx*sx; a=(sy*sx2-sx*sxy)/deno; b=(n*sxy-sy*sx)/deno; printf(a=%6.2f b=%6.2fn,a,b); return(0);= 例题3.2.2计算

10、程序! nonlinearfit.for program nonlinearfit! non-linear_fitimplicit none integer n,i,mparameter (n=8,m=20)real, dimension (n): x,y,yc real a,b,dxdata x/1,2,3,4,5,6,7,8/data y/15.3,20.5,27.4,36.6,49.1,65.6,87.8,117.6/ do i=1,n yc(i)=log(y(i)end docall fit(n,x,yc,a,b)a=exp(a)print *,a=,a, b=,b endsubrou

11、tine fit(n,x,y,a,b)implicit noneinteger n,ireal, dimension (n): x,y real xp,yp,x2,xyp,a,b,cxp = 0.0yp = 0.0x2 = 0.0xyp= 0.0do i = 1,n x2 = x2+x(i)*x(i) xp = xp+x(i) yp = yp+y(i) xyp= xyp+x(i)*y(i)end do xp = xp/nyp = yp/nx2 = x2/nxyp= xyp/nc = x2-xp*xp a = (yp*x2-xp*xyp)/cb = (xyp-xp*yp)/c return en

12、d=3.2.3 m次多项式拟合 例题3.2.3计算程序! polyfit.For program polyfit! 多项式拟合implicit noneinteger n,m,i,l,k,jparameter (n=11,m=6,k=21) real, dimension (n): x,yreal, dimension (k): xn,yn,xo,yo real s(m,m),t(m),js(m),z(m),dxopen(5,file=exa3_2_3.dat) open(2,file=exa3_2_3_old.dat) data x/-1.0,-0.8,-0.6,-0.4,-0.2,0.0,

13、0.2,0.4,0.6,0.8,1.0/do i=1,n y(i)=exp(x(i)end docall fitp(n,m,x,y,s,t) call gaus(s,t,m,z,l,js)if (l.ne.0) then write(*,10) (i,z(i),i=1,6)end if10format(1x,x(,i2, )=,d15.6) dx=2.0/(k-1)do i=1,k xn(i)=-1.0+dx*(i-1) xo(i)=xn(i) yo(i)=exp(xo(i) yn(i)=z(m) do j=1,m-1 yn(i)=yn(i)*xn(i)+z(m-j) end do write

14、(5,(2f14.6) xn(i),yn(i) write(2,(2f14.6) xo(i),yo(i)end doendsubroutine fitp(n,m1,x,y,s,t)implicit noneinteger n,m1,m,j,i,i1,mireal, dimension (n): x,y real s(m1,m1),t(m1)m=m1-1s=0.0t=0.0s(1,1)=ndo j=1,n t(1)=t(1)+y(j)end dodo i=2,m1 i1=i-1 mi=m1+i-2 do j=1,n s(i,1)=s(i,1)+x(j)*i1 s(m1,i)=s(m1,i)+x(

15、j)*mi t(i)=t(i)+x(j)*i1*y(j) end doend dodo j=2,m do i=j,m s(i,j)=s(i+1,j-1) end doend dodo i=1,m i1=i+1 do j=i1,m1 s(i,j)=s(j,i) end doend doreturnend subroutine gaus(a,b,n,x,l,js)real a(n,n),x(n),b(n),js(n)real tl=1do 50 k=1,n-1 d=0.0 do 210 i=k,n do 210 j=k,n if (abs(a(i,j).gt.d) then d=abs(a(i,j

16、) js(k)=j is=i end if210 continue if (d+1.0.eq.1.0) then l=0 else if (js(k).ne.k) then do 220 i=1,n t=a(i,k) a(i,k)=a(i,js(k) a(i,js(k)=t220 continue end if if (is.ne.k) then do 230 j=k,n t=a(k,j) a(k,j)=a(is,j) a(is,j)=t230 continue t=b(k) b(k)=b(is) b(is)=t end if end if if (l.eq.0) then write(*,1

17、00) return end if do 10 j=k+1,n a(k,j)=a(k,j)/a(k,k)10 continue b(k)=b(k)/a(k,k) do 30 i=k+1,n do 20 j=k+1,n a(i,j)=a(i,j)-a(i,k)*a(k,j)20 continue b(i)=b(i)-a(i,k)*b(k)30 continue50continueif (abs(a(n,n)+1.0.eq.1.0) then l=0 write(*,100) returnend ifx(n)=b(n)/a(n,n)do 70 i=n-1,1,-1 t=0.0 do 60 j=i+

18、1,n t=t+a(i,j)*x(j)60 continue x(i)=b(i)-t70continue100format(1x, fail )js(n)=ndo 150 k=n,1,-1 if (js(k).ne.k) then t=x(k) x(k)=x(js(k) x(js(k)=t end if150continuereturnend% demo_polyfitsclc;clear all;format long;x =-1:0.2:1;y = exp(x);m = 5;a = polyfits(x,y,m)c =a(6) a(5) a(4) a(3) a(2) a(1);x1=-1:0.02:1;y1=polyval(c,x1);plot(x1,y1,r,x,y,o);title(5次多项式拟合 ex