《数值分析插值法.docx》由会员分享，可在线阅读，更多相关《数值分析插值法.docx（14页珍藏版）》请在三一文库上搜索。
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)=Yfor j=2:nfor k=j:nD(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:1C=conv(C,p
2、oly(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;X=-1:0.1:1;Y=1./(1+25*
3、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(2:N);U=6*diff(D);B(1)=B(1)-H(1)/2
4、;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-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);for k=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;for k=0:N-1S(k+1,1)=(M(k+2)-M(k
5、+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,:),x
6、2-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)
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=polyf
9、it(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);结果如图：（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+1if k=jV=conv(V,poly(x(j)/(x(k)-x(j);endendL(k,:)=V;endC=y*L在命令窗口中输入以下的命令：x=0.
10、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);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.0
11、1:1.0;yy=polyval(C,xx);hold onplot(xx,yy,b,x,y,.);第五章解线性方程组的直接方法1.用 LU 分解及列主元消去法解线性方程组10701x183 2.09999962x25.9000015151x352103x41输出 Ax=b 中系数 A=LU 分解的矩阵L 及 U，解向量x 及 detA；列主元法的行交换次序，解向量 x 及 detA；比较两种方法所得的结果。解：程序如下：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)输
12、出的结果如下：求 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 = 1;for j = 1: (m - 1)if (A(i, j) = 0)for k = (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, :) =
13、 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);for i = (n - 1):-1:1temp_b = 0;for j = (i + 1):ntemp_b = temp_b + A(i, j) * X(j);endX(i) = (b(i) - temp_b) / A(i, i);endend输出结果为：A=10 -7 0 1;-3 2.099999 6 2;
14、5 -1 5 -1;2 1 0 2第八章矩阵特征值的计算1078775652345611/ 21/ 31/ 41/ 51/ 686109445671/ 21/ 31/ 41/ 51/ 61/ 71.已知矩阵 A= 75910 ,B= 03678,= 1/ 31/ 41/ 51/ 61/ 71/ 8002891/ 41/ 51/ 61/ 71/ 81/ 9000101/ 51/ 61/ 71/ 81/ 91/101/ 61/ 71/ 81/ 91/101/11（ 1）用 MATLAB 函数“ eig”求矩阵全部特征值。（ 2）用基本 QR 算法求全部特征值（可用 MATLAB函数“ qr”实现矩
15、阵的 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 的特征值：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)输出结果：求矩阵的特征值：clear;=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/81/9;1/5 1/6 1/7 1/8
16、 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11;E=eig()输出结果：10787（2） A= 75658610975910第一步： A0=hess(A);Q0,R0=qr(A0);A1=R0*Q0返回得到：第二部： Q1,R1=qr(A1);A2=R1*Q1第三部： Q2,R2=qr(A2);A3=R2*Q229.83293.44110.0000现在收缩，继续对A3 的子矩阵=3.44114.30530.161100.16110.8516进行累世变换，得到（假设收缩后的矩阵为C6）29.83293.44110.0000C6=3.44114.30530.161100.16110.8516这是进行了6 步 qr 算法所得的结果。故求的A 的近似特征值为30.2886 ，0.0102。而 A 的特征值是0.010230.2886同理，用类似的方法可求矩阵B 和的特征值，但过程过于繁琐，不再一一求解。