(完整版)数值分析插值法

第二章插值法
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:n
for k=j:n
D(k,j)=(D(k,j-1)- D(k-1,j-1))/(X(k)-X(k-j+1));
end
end
C=D(n,n);
for k=(n-1):-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;
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-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);
end
M(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-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,clc
X=-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,clc
X=-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. 由实验给出数据表
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
试求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+1
V=1;
for j=1:n+1
if k~=j
V=conv(V,poly(x(j)))/(x(k)-x(j));
end
end
L(k,:)=V;
end
C=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 on
plot(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):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;
A
b
disp('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: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);
end
end
输出结果为：
A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]
第八章矩阵特征值的计算
1.已知矩阵A=
10787
7565
86109
75910
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦,B=
23456
44567
03678
00289
00010
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
,=
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/51/61/71/81/91/10
1/61/71/81/91/101/11
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎢⎥
⎣⎦
（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的特征值：
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/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()
输出结果：
（2）A=10787 7565 86109 75910
第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到：
第二部：[Q1,R1]=qr(A1);A2=R1*Q1
第三部：[Q2,R2]=qr(A2);A3=R2*Q2
现在收缩，继续对A3的子矩阵=29.8329 3.44110.0000 3.4411 4.30530.1611
00.16110.8516
-
-
进行累世变换，得到（假设收缩后的矩阵为C6）
C6=29.8329 3.44110.0000 3.4411 4.30530.1611
00.16110.8516
-
-
这是进行了6步qr算法所得的结果。

故求的A的近似特征值为30.2886，，，0.0102。

而A的特征值是0.0102
30.2886
同理，用类似的方法可求矩阵B 和的特征值，但过程过于繁琐，不再一一求解。