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

1第3章附录3.1.5n+1点n次插值functiony=lagrange_2(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=s+p*y0(k);endy(i)=s;End==============================================3.1.6三次样条插值例题3.1.7计算程序%demo_spline.mx=-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.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.)fori=2:n-1a(i)=x(i+1)-x(i);b(i)=2.*(x(i+1)-x(i));2c(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)));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);fori=1:n-1c1(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));endj=1;fori=1:n-1dj=floor((x(i+1)-x(i)+dx/2)/dx);fork=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);endj=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('原曲线','插值点','插值曲线');title('y=1/(1+x^2):样条插值','FontSize',12);==============================================例题3.2.1计算程序!!!!linear_fit.for!!!programlinear_fit!linear_fitimplicitnoneintegern,i,mparameter(n=7,m=20)real,dimension(n)::x,yreal,dimension(m)::xi,yireala,b,dx3open(5,file='exa3_2_1.dat')open(2,file='exa3_2_1_old.dat')datax/0.5,1.2,2.1,2.9,3.6,4.5,5.7/datay/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)callfit(n,x,y,a,b)print*,'a=',a,'b=',bdx=(x(n)-x(1))/(m-1)doi=1,mxi(i)=x(1)+(i-1)*dxyi(i)=a+b*xi(i)enddowrite(2,'(2f14.6)')(xi(i),yi(i),i=1,m)endsubroutinefit(n,x,y,a,b)implicitnoneintegern,ireal,dimension(n)::x,yrealxp,yp,x2,xyp,a,b,cxp=0.0yp=0.0x2=0.0xyp=0.0doi=1,nx2=x2+x(i)*x(i)xp=xp+x(i)yp=yp+y(i)xyp=xyp+x(i)*y(i)enddoxp=xp/nyp=yp/nx2=x2/nxyp=xyp/nc=x2-xp*xpa=(yp*x2-xp*xyp)/cb=(xyp-xp*yp)/creturnend4//*linearfit.cpp*/#includestdio.hmain(){intn=7,i;floatx[]={0.5,1.2,2.1,2.9,3.6,4.5,5.7};floaty[]={2.81,3.24,3.80,4.30,4.73,5.29,6.03};floatsx=0.,sy=0.,sxy=0.,sx2=0,deno,a,b;for(i=0;i=6;i++){sx=sx+x[i];sy=sy+y[i];sxy=sxy+x[i]*y[i];sx2=sx2+x[i]*x[i];}deno=n*sx2-sx*sx;a=(sy*sx2-sx*sxy)/deno;b=(n*sxy-sy*sx)/deno;printf(a=%6.2fb=%6.2f\n,a,b);return(0);}==============================================例题3.2.2计算程序!nonlinearfit.forprogramnonlinearfit!non-linear_fitimplicitnoneintegern,i,mparameter(n=8,m=20)real,dimension(n)::x,y,ycreala,b,dxdatax/1,2,3,4,5,6,7,8/datay/15.3,20.5,27.4,36.6,49.1,65.6,87.8,117.6/doi=1,nyc(i)=log(y(i))enddocallfit(n,x,yc,a,b)a=exp(a)print*,'a=',a,'b=',bend5subroutinefit(n,x,y,a,b)implicitnoneintegern,ireal,dimension(n)::x,yrealxp,yp,x2,xyp,a,b,cxp=0.0yp=0.0x2=0.0xyp=0.0doi=1,nx2=x2+x(i)*x(i)xp=xp+x(i)yp=yp+y(i)xyp=xyp+x(i)*y(i)enddoxp=xp/nyp=yp/nx2=x2/nxyp=xyp/nc=x2-xp*xpa=(yp*x2-xp*xyp)/cb=(xyp-xp*yp)/creturnend=============================================3.2.3m次多项式拟合例题3.2.3计算程序!!!polyfit.Forprogrampolyfit!多项式拟合implicitnoneintegern,m,i,l,k,jparameter(n=11,m=6,k=21)real,dimension(n)::x,yreal,dimension(k)::xn,yn,xo,yoreals(m,m),t(m),js(m),z(m),dxopen(5,file='exa3_2_3.dat')open(2,file='exa3_2_3_old.dat')datax/-1.0,-0.8,-0.6,-0.4,-0.2,0.0,0.2,0.4,0.6,0.8,1.0/doi=1,ny(i)=exp(x(i))enddo6callfitp(n,m,x,y,s,t)callgaus(s,t,m,z,l,js)if(l.ne.0)thenwrite(*,10)(i,z(i),i=1,6)endif10format(1x,'x(',i2,')=',d15.6)dx=2.0/(k-1)doi=1,kxn(i)=-1.0+dx*(i-1)xo(i)=xn(i)yo(i)=exp(xo(i))yn(i)=z(m)doj=1,m-1yn(i)=yn(i)*xn(i)+z(m-j)enddowrite(5,'(2f14.6)')xn(i),yn(i)write(2,'(2f14.6)')xo(i),yo(i)enddoendsubroutinefitp(n,m1,x,y,s,t)implicitnoneintegern,m1,m,j,i,i1,mireal,dimension(n)::x,yreals(m1,m1),t(m1)m=m1-1s=0.0t=0.0s(1,1)=ndoj=1,nt(1)=t(1)+y(j)enddodoi=2,m1i1=i-1mi=m1+i-2doj=1,ns(i,1)=s(i,1)+x(j)**i1s(m1,i)=s(m1,i)+x(j)**mit(i)=t(i)+x(j)**i1*y(j)enddoenddo7doj=2,mdoi=j,ms(i,j)=s(i+1,j-1)enddoenddodoi=1,mi1=i+1doj=i1,m1s(i,j)=s(j,i)enddoenddoreturnendsubroutinegaus(a,b,n,x,l,js)reala(n,n),x(n),b(n),js(n)realtl=1do50k=1,n-1d=0.0do210i=k,ndo210j=k,nif(abs(a(i,j)).gt.d)thend=abs(a(i,j))js(k)=jis=iendif210continueif(d+1.0.eq.1.0)thenl=0elseif(js(k).ne.k)thendo220i=1,nt=a(i,k)a(i,k)=a(i,js(k))a(i,js(k))=t220continueendifif(is.ne.k)thendo230j=k,nt=a(k,j)a(k,j)=a(is,j)8a(is,j)=t230continuet=b(k)b(k)=b(is)b(is)=tendifendifif(l.eq.0)thenwrite(*,100)returnendifdo10j=k+1,na(k,j)=a(k,j)/a(k,k)10continueb(k)=b(k)/a(k,k)do30i=k+1,ndo20j=k+1,na(i,j)=a(i,j)-a(i,k)*a(k,j)20continueb(i)=b(i)-a(i,k)*b(k)30continue50continueif(abs(a(n,n))+1.0.eq.1.0)thenl=0write(*,100)returnendifx(n)=b(n)/a(n,n)do70i=n-1,1,-1t=0.0do60j=i+1,nt=t+a(i,j)*x(j)60continuex(i)=b(i)-t70co