数值分析-插值法

第二章插值法2.在区间[-1,1]上分别取n=10,20用两组等距节点对龙哥函数f(x)=1/(1+25*x^2)做多项式插值及三次样条插值，对每个n值，分别画出插值函数及f（x）的图形。（1）多项式插值①先建立一个多项式插值的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.^2);[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.^2);plot(xp,z,'r')得到插值函数和f（x）图形：③当n=20时，我们在命令窗口中输入以下的命令：clear,clf,holdon;X=-1:0.1:1;Y=1./(1+25*X.^2);[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.^2);plot(xp,z,'r')得到插值函数和f（x）图形：（2）三次样条插值①先建立一个多项式插值的M-file；输入如下的命令:functionS=csfit(X,Y,dx0,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.^2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Y,dx0,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(2));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.^2+1);dx0=0.0739644970414201;dxn=-0.0739644970414201;S=csfit(X,Y,dx0,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(2));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');结果如图：（3）、另一个拟合曲线：新建一个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分解及列主元消去法解线性方程组123410701832.099999625.9000015151521031xxxx输出Ax=b中系数A=LU分解的矩阵L及U，解向量x及detA；列主元法的行交换次序，解向量x及detA；比较两种方法所得的结果。解：程序如下：clear;A=[10-701;-32.09999962;5-15-1;2102];B=[8;5.900001;5;1];[L,U]=lu(A);X=U\(L\B)输出的结果如下：求det（A）：输入：det(A);输出：列主元素消去法：程序如下：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]第八章矩阵特征值的计算1.已知矩阵A=1078775658610975910,B=2345644567036780028900010,=11/21/31/41/51/61/21/31/41/51/61/71/31/41/51/61/71/81/41/51/61/71/81/91/51/61/71/81/91/101/61/71/81/91/101/11（1）用MATLAB函数“eig”求矩阵全部特征值。（2）用基本QR算法求全部特征值（可用MATLAB函数“qr”实现矩阵的QR分解）。解：MATLAB程序如下：求矩阵A的特征值：clear;A=[10787;7565;86109;75910];E=eig(A)输出结果：求矩阵B的特征值：clear;B=[23456;44567;03678;00289;00010];E=eig(B)输出结果：求矩阵的特征值：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/51/61/71/81/91/10;1/61/71/81/91/101/11];E=eig()输出结果：（2）A=1078775658610975910第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0返回得到：第二部：[Q1,R1]=qr(A1);A2=R1*Q1第三部：[Q2,R2]=qr(A2);A3=R2*Q2现在收缩，继续对A3的子矩阵=29.83293.44110.00003.44114.30530.161100.16110.8516进行累世变换，得到（假设收缩后的矩阵为C6）C6=29.83293.44110.00003.44114.30530.161100.16110.8516这是进行了6步qr算法所得的结果。故求的A的近似特征值为30.2886，，，0.0102。而A的特征值是0.010230.2886同理，用类似的方法可求矩阵B和的特征值，但过程过于繁琐，不再一一求解。