第二章插值法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