數(shù)值分析-插值法

1、第二章 插值法2.在區(qū)間-1,1上分別取n=10,20用兩組等距節(jié)點(diǎn)對(duì)龍哥函數(shù) f(x)=1/(1+25*x2)做多項(xiàng)式插值及三次樣條插值，對(duì)每個(gè)n值，分別畫(huà)出插值函數(shù)及f（x）的圖形。（1）多項(xiàng)式插值先建立一個(gè)多項(xiàng)式插值的M-file；輸入如下的命令（如牛頓插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y' for j=2:n for k=j:n D(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1); endendC=D(n,n);for k=(n-1):-1:1 C=c

2、onv(C,poly(X(k) m=length(C); C(m)= C(m)+D(k,k);end當(dāng)n=10時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clf,hold on; 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,'.'); grid on; xp=-1:0.2:1; z=1./(1+25*xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：當(dāng)n=20時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,

3、clf,hold on; 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,'.'); grid on; xp=-1:0.1:1; z=1./(1+25*xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：（2）三次樣條插值先建立一個(gè)多項(xiàng)式插值的M-file；輸入如下的命令:function S=csfit(X,Y,dx0,dxn) N=length(X)-1;H=diff(X); D=diff(Y)./H;A=H(2

4、: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); for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1);for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1

5、)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0:N-1 S(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當(dāng)n=10時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1

6、=-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,'.')結(jié)果如圖：當(dāng)n=20時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clcX=-1:0.1:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;d

7、xn= -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,'.')結(jié)果如圖：第三章 函數(shù)逼近與快速傅里葉變換2. 由實(shí)驗(yàn)給出數(shù)據(jù)表x0.00.10.20

8、.30.50.81.0y1.00.410.500.610.912.022.46試求3次、4次多項(xiàng)式的曲線擬合，再根據(jù)數(shù)據(jù)曲線形狀，求一個(gè)另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。（1）、三次擬合曲線：命令如下：x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x,y,3);xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,'-');hold on;plot(x,y,'x');xlabel(&

9、#39;x');ylabel('y');結(jié)果如圖：（2）、4次擬合曲線輸入命令：x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,'r');hold on;plot(x,y,'x');xlabel('x');ylabel('y');結(jié)果如圖：（3）、另一個(gè)擬合曲線：新建一個(gè)M-file：

10、輸入如下命令：function C,L=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);for k=1:n+1V=1;for j=1:n+1 if k=jV=conv(V,poly(x(j)/(x(k)-x(j);endend L(k,:)=V;endC=y*L在命令窗口中輸入以下的命令：x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=1.0 0.41 0.50 0.61 0.91 2.02 2.46;cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,

11、'r');hold on;plot(x,y,'x');xlabel('x');ylabel('y');x=0.0 0.1 0.2 0.3 0.5 0.8 1.0;y=0.94 0.58 0.47 0.52 1.00 2.00 2.46; %y中的值是根據(jù)上面兩種擬合曲線而得到的估計(jì)數(shù)據(jù)，不是真實(shí)數(shù)據(jù)C,L=lagran(x,y); xx=0:0.01:1.0;yy=polyval(C,xx);hold onplot(xx,yy,'b',x,y,'.');第五章 解線性方程組的直接方法1.用LU分解及

12、列主元消去法解線性方程組輸出Ax=b中系數(shù)A=LU分解的矩陣L及U，解向量x及detA；列主元法的行交換次序，解向量x及detA；比較兩種方法所得的結(jié)果。解：程序如下：clear;A=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2;B=8;5.900001;5;1;L,U=lu(A);X=U(LB)輸出的結(jié)果如下：求det（A）：輸入：det(A);輸出：列主元素消去法：程序如下：function X =Gauss(A, b)n, m = size(A);X = zeros(n, 1);temp = zeros(1, m);temp_b = 0;i =

13、1;for j = 1: (m - 1) if (A(i, j) = 0) for k = (i + 1):n if (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; end end end i = i + 1;end;Abdisp('det(A) is .');x = det(A);disp(x);disp('cond(A) is .');

14、x = cond(A);disp(x);X(n) = b(n) / A(n, n);for i = (n - 1):-1:1 temp_b = 0; for j = (i + 1):n temp_b = temp_b + A(i, j) * X(j); end X(i) = (b(i) - temp_b) / A(i, i);endend輸出結(jié)果為：A=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2第八章 矩陣特征值的計(jì)算1.已知矩陣A=,B=,H6=（1）用MATLAB函數(shù)“eig”求矩陣全部特征值。（2）用基本QR算法求全部特征值（可用MATLAB函

15、數(shù)“qr”實(shí)現(xiàn)矩陣的QR分解）。解：MATLAB程序如下：求矩陣A的特征值：clear;A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;E=eig(A)輸出結(jié)果：求矩陣B的特征值：clear;B=2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0;E=eig(B)輸出結(jié)果：求矩陣H6的特征值：clear;H6=1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11;E=eig(H6)輸出結(jié)果：