《数值分析-插值法(共13页)》由会员分享，可在线阅读，更多相关《数值分析-插值法(共13页)（13页珍藏版）》请在金锄头文库上搜索。
1、精选优质文档-倾情为你奉上第二章 插值法2.在区间-1,1上分别取n=10,20用两组等距节点对龙哥函数 f(x)=1/(1+25*x2)做多项式插值及三次样条插值，对每个n值，分别画出插值函数及f（x）的图形。（1）多项式插值先建立一个多项式插值的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):
2、-1:1 C=conv(C,poly(X(k) m=length(C); C(m)= C(m)+D(k,k);end当n=10时，我们在命令窗口中输入以下的命令：clear,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)得到插值函数和f（x）图形：当n=20时，我们在命令窗口中输入以下的命令：clear,clf,hold on;
3、 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)得到插值函数和f（x）图形：（2）三次样条插值先建立一个多项式插值的M-file；输入如下的命令:function S=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(
4、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)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)
5、-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当n=10时，我们在命令窗口中输入以下的命令：clear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.14201;dxn= -0.14201;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
6、;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.14201;dxn= -0.14201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1)
7、;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次多项式的曲线拟合，再根据数据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。
8、（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(x);ylabel(y);结果如图：（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);
9、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);结果如图：（3）、另一个拟合曲线：新建一个M-file：输入如下命令：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
10、 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);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中的值是根据上面两种拟合曲线而得到的估计数据，不是真实数据C,L=lagran(x,y); xx=0:0.01:1.0;yy=
11、polyval(C,xx);hold onplot(xx,yy,b,x,y,.);第五章 解线性方程组的直接方法1.用LU分解及列主元消去法解线性方程组输出Ax=b中系数A=LU分解的矩阵L及U，解向量x及detA；列主元法的行交换次序，解向量x及detA；比较两种方法所得的结果。解：程序如下：clear;A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;B=8;5.;5;1;L,U=lu(A);X=U(LB)输出的结果如下：求det（A）：输入：det(A);输出：列主元素消去法：程序如下：function X =Gauss(A, b)n, m = size(A
12、);X = zeros(n, 1);temp = zeros(1, m);temp_b = 0;i = 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 = de
13、t(A);disp(x);disp(cond(A) is .);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输出结果为：A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2第八章 矩阵特征值的计算1.已知矩阵A=,B=,=（1）用MATLAB函数“eig”求矩阵全部特征值。（2）用基本QR算法求全部特征值（可用MATLAB函数“qr”实现矩阵的QR分解）。解：MATLAB程序如下：求矩阵A的特征值：clear;A=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;E=eig(A)输出结果：求矩阵B的特征值：c
