数值分析 插值法

第二章插值法2.在区间[-1,1]上分别取n=10,20用两组等距节点对龙哥函数f(x)=1/(1+25*xA2)做多项式插值及三次样条插值，对每个n值，分别画出插值函数及f(x)的图形。多项式插值先建立一个多项式插值的M-file；输入如下的命令(如牛顿插值公式)：function[C,D]=newpoly(X,Y)n=length(X);D=zeros(n,n)D(:,1)=Y'forj=2:nfork=j:nD(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);fork=(n-1):-1:1C=conv(C,poly(X(k)))m=length(C);C(m)=C(m)+D(k,k);end当n=10时，我们在命令窗口中输入以下的命令：clear,clf,holdon;X=-1:0.2:1;Y=1./(1+25*X.A2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y'.');gridon;xp=-1:0.2:1;z=1./(1+25*xp.A2);plot(xp,z,'r')得到插值函数和f(x)图形：

当n=20时，我们在命令窗口中输入以下的命令:clear,clf,holdon;X=-1:0.1:1;Y=1./(1+25*X.A2);[C,D]=newpoly(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y'.');gridon;xp=-1:0.1:1;z=1./(1+25*xp.A2);plot(xp,z,'r')得到插值函数和f(x)图形：三次样条插值先建立一个多项式插值的M-file；输入如下的命令：functionS=csfit(X,Ydx0,dxn)N=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N));C=H(2:N);U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1));B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N));fork=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);fork=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;fork=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);end当n=10时，我们在命令窗口中输入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.A2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Ydx0,dxn)x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X⑵)；x3=0:0.01:0.5;y3=polyval(S(3,:),x3-X(3));x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4));plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y'.')结果如图：

②当n=20时，我们在命令窗口中输入以下的命令：clear,clcX=-1:0.1:1;Y=1./(25*X.A2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Ydx0,dxn)x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1));x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X⑵)；x3=0:0.01:0.5;y3=polyval(S(3,:),x3-X(3));x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4));plot(x1,y1,x2,y2,x3,y3,x4,y4,X,Y'.')结果如图：

第三章函数逼近与快速傅里叶变换2.由实验给出数据表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46试求3次、4次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。(1)、三次拟合曲线：命令如下：x=[0.00.10.20.30.50.81.0];y=[1.00.410.500.610.912.022.46];cc=polyfit(x,y,3);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'--');holdon;plot(x,y,'x');xlabel('x');ylabel('y');结果如图：(2)、4次拟合曲线输入命令：x=[0.00.10.20.30.50.81.0];y=[1.00.410.500.610.912.022.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');holdon;plot(x,y,'x');xlabel('x');ylabel('y');结果如图：、另一个拟合曲线：新建一个M-file：输入如下命令：function[C,L]=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);fork=1:n+1V=1;forj=1:n+1ifk~=jV=conv(V,poly(x(j)))/(x(k)-x(j));endendL(k,:)=V;endC=y*L在命令窗口中输入以下的命令：x=[0.00.10.20.30.50.81.0];y=[1.00.410.500.610.912.022.46];cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x));yy=polyval(cc,xx);plot(xx,yy,'r');holdon;plot(x,y,'x');xlabel('x');ylabel('y');x=[0.00.10.20.30.50.81.0];y=[0.940.580.470.521.002.002.46];%y中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据[C,L]=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdonplot(xx,yy,'b',x,y,'.');第五章解线性方程组的直接方法1.用LU分解及列主元消去法解线性方程组10-352-72.099999061-2「尤]1X2—_「8一5.900001-110-352-72.099999061-2「尤]1X2—_「8一5.900001-15-1X35103X1-4」1解：程序如下：clear;A=[10-701;-32.09999962;5-15-1;2102];B=[8;5.900001;5;1];[L,U]=lu(A);X=U\(L\B)输出的结果如下：>>clearA=[10-701;-32.09999962;5-15-1;2102];[L,U]=lu(A)[L,U]=lu(A)000-0.3000-0.00001.000000.50001.0000000.20000.9600-0.80001.0000u=10.0000-7.000001.000002.50005.0000-1.5000006.00002.30000005.0800»B=[8;5.900001;5;1];X=U\(L\B)0.0000-1.00001.00001.0000求det（A）：输入：det(A);输出：»det(A)ans=-762.0001列主元素消去法：程序如下：functionX=Gauss(A,b)[n,m]=size(A);X=zeros(n,1);temp=zeros(1,m);temp_b=0;i=1;forj=1:(m-1)if(A(i,j)~=0)fork=(i+1):nif(A(k,j)~=0)temp=A(k,:)+A(i,:)*(-A(k,j)/A(i,j));temp_b=b(k)+b(i)*(-A(k,j)/A(i,j));A(k,:)=temp;b(k)=temp_b;endendendi=i+1;end;Abdisp('det(A)is...');x=det(A);disp(x);disp('cond(A)is...');x=cond(A);disp(x);X(n)=b(n)/A(n,n);fori=(n-1):-1:1temp_b=0;forj=(i+1):ntemp_b=temp_b+A(i,j)*X(j);endX(i)=(b(i)-temp_b)/A(i,i);endend输出结果为：A=[10-701;-32.09999962;5-15-1;2102]10.0000-7.000001.0000-3.00002.10006.00002.00005.0000-1.00005.0000-1.00002.00001.000002.0000»d=[8;5.900001;5;1]d二8.00005.90005.00001.0000>>z=Gauss(Ajd)k=-0.0000-1.00001.00001.0000第八章矩阵特征值的计算1.已知矩阵A=1078775658610975910,B=6一一11/21/31/41/51/6-71/21/31/41/51/61/78丹=1/31/41/51/61/71/891/41/51/61/71/81/901/51/61/71/81/91/10—_1/61/71/81/91/101/11用MATLAB函数“eig”求矩阵全部特征值。用基本QR算法求全部特征值(可用MATLAB函数“qr”实现矩阵的QR分解)。解：MATLAB程序如下：求矩阵A的特征值：clear;A=[10787;7565;86109;75910];E=eig(A)输出结果：»A=[10787;7565;86109;75910];E=eig(A)0.01020.84313.858130.2887求矩阵B的特征值：clear;B=[23456;44567;03678;00289;00010];E=eig(B)输出结果：>>clear»B=[23456;44567;03678;00289;00010];E=eig(B)E二13.17246.55191.5957-0.3908-0.9291求矩阵％的特征值:clear;=[11/21/31/41/51/6;1/21/31/41/51/61/7;1/31/41/51/61/71/8;1/41/51/61/71/81/9;1/5161/71/81/91/10;161/71/81/91/101/11];E=eig(H(?)输出结果：>>clear»H6=[l1/21/31/41/51/6;1/21/31/41/51/61/7;1/31/41/51/61/71/8;1/41/51/61/71/81/9;1/51/61/71/81/91/10;1/61/71/81/91/101/11]:E=eig(H6)E二0.l：ll：ll：ll：l0.i：ii：ii：ii：i0.0i：ii：i60.01630.24241.6189TOC\o"1-5"\h\z10787(2)A=75658610975910第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0返回得到：A1二2.5597-2.83640.0000-0.0000-2.836422.178212.51870.0000012.518710.24860.0808000.08080.0134第二部：[Q1,R1]=qr(A1);A2=R1*Q1A2=16.1940-13.26690.00000.0000-13.266917.76611.0092-0.000001.00921.0297-0.000600-0.00060.0102第三部：[Q2,R2]=qr(A2);A3=R2*Q2心=29.8329-3.44110.0000-0.0000-3.44114.30530.1611-0