当前位置:文档之家› 数值分析 插值法

数值分析 插值法

第二章插值法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'for 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,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;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(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-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+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.试求3次、4次多项式的曲线拟合,再根据数据曲线形状,求一个另外函数的拟合曲线,用图示数据曲线及相应的三种拟合曲线。

(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);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.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.01:1.0;yy=polyval(C,xx);hold onplot(xx,yy,'b',x,y,'.');第五章 解线性方程组的直接方法1.用LU 分解及列主元消去法解线性方程组12341070183 2.09999962 5.9000015151521031x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦输出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\(L\B) 输出的结果如下:求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, :) = 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;5 -1 5 -1;2 1 0 2]第八章矩阵特征值的计算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”求矩阵全部特征值。

相关主题