矩阵分析与数值分析实验报告

《矩阵分析与数值分析》实验报告
院系：
姓名：
学号：
所在班号：
任课老师：
一．设错误！未找到引用源。

，分别编制从小到大和从大到小的顺序程序计算错误！未找到引用源。

并指出有效位数。

程序如下：
function sum3
j=input('请输入求和个数 "j":');
A=0;
B=0;
double B;
double A;
for n=2:j
m=n^2-1;
t=1./m;
A=A+t;
end
disp('从小到大：')
s=A
for n=j:-1:2
m=n^2-1;
t=1./m;
B=B+t;
end
disp('从大到小:')
s=B
运行结果：
>> sum3
请输入求和个数 "j":100
从小到大：
s =0.740049504950495
从大到小:
s =0.740049504950495
>> sum3
请输入求和个数 "j":10000
从小到大：
s =0.749900004999506
从大到小:
s =0.749900004999500
>> sum3
请输入求和个数 "j":1000000
从小到大：
s =0.749999000000522
从大到小:
s =0.749999000000500
二、解线性方程组
1．分别Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组。

⎪⎪⎪⎪⎪⎭
⎫
⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪
⎪
⎪⎪
⎭
⎫ ⎝
⎛----000121001210
0121
00124321x x x x 迭代法计算停止的条件为：6)()1(3
110max -+≤≤<-k j k j j x x 。

解：（1）Jacobi 迭代法程序代码： function jacobi(A, b, N) clc;clear;
A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2]; b=[-1 0 0 0]'; N=100;
n = size(A,1); D = diag(diag(A)); L = tril(-A,-1); U = triu(-A,1); Tj = inv(D)*(L+U); cj = inv(D)*b; tol = 1e-06; k = 1;
format long
x = zeros(n,1); while k <= N
x(:,k+1) = Tj*x(:,k) + cj;
disp(k); disp('x = ');disp(x(:,k+1)); if norm(x(:,k+1)-x(:,k)) < tol
disp('The procedure was successful')
disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1)); break end
k = k+1; end
结果输出
The procedure was successful
Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations 60 x =
0.79999879906731
0.59999842795870
0.39999805685009
0.19999902842505
（2）Gauss-Seidel迭代法程序代码：
function gauss_seidel(A, b, N)
clc;clear;
A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];
b=[-1 0 0 0]';
N=100;
n = size(A,1);
D = diag(diag(A));
L = tril(-A,-1);
U = triu(-A,1);
Tg = inv(D-L)*U;
cg = inv(D-L)*b;
tol = 1e-06;
k = 1;
x = zeros(n,1);
while k <= N
x(:,k+1) = Tg*x(:,k) + cg;
disp(k); disp('x = ');disp(x(:,k+1));
if norm(x(:,k+1)-x(:,k)) < tol
disp('The procedure was successful')
disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1));
break
end
k = k+1;
end
结果输出
The procedure was successful
Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations
31
x =
0.79999921397935
0.59999897108561
0.39999916759077
0.19999958379539
2. 用Gauss列主元消去法、QR方法求解如下方程组：
⎪⎪⎪⎪⎪⎭
⎫
⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪
⎪
⎪⎪
⎭
⎫ ⎝⎛---017435222331
325212214321x x x x (1)Gauss 列主元消去法 程序代码：
function x=Gaussmain(A,b) clc;clear; format long
A=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3]; b=[4 7 -1 0]'; N=length(A); x=zeros(N,1); y=zeros(N,1); c=0; d=0;
A(:,N+1)=b; for k=1:N-1 for i=k:4
if c<abs(A(i,k))
d=i;c=abs(A(i,k)); end end
y=A(k,:);
A(k,:)=A(d,:); A(d,:)=y; for i=k+1:N c=A(i,k);
for j=1:N+1
A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end end end
b=A(:,N+1);
x(N)=b(N)/A(N,N); for k=N-1:-1:1
x(k)=(b(k)-A(k,k+1:N)*x(k+1:N))/A(k,k); end
结果输出 ans =
18.00000000000000 -9.57142857142857
6.00000000000000
-0.42857142857143
（2）QR方法：程序代码
function QR(A,b)
clc;clear;
format long
A=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3];
b=[4 7 -1 0]';
[Q,R]=qr(A)
y=Q'*b;
x=R\y
结果输出
Q =
-0.31622776601684 -0.04116934847963 -0.75164602800283 0.57735026918962 -0.63245553203368 -0.49403218175557 -0.15032920560056 -0.57735026918963 0.63245553203368 -0.74104827263336 -0.22549380840085 -0.00000000000000 -0.31622776601684 -0.45286283327594 0.60131682240226 0.57735026918963 R =
-3.16227766016838 -6.00832755431992 -0.94868329805051 2.84604989415154 0 -2.42899156029822 -4.65213637819829 -4.15810419644272 0 0 -0.67648142520255 -0.52615221960200 0 0 0 4.04145188432738 x =
17.99999999999989
-9.57142857142851
5.99999999999997
-0.42857142857143
三、非线性方程的迭代解法
1．用Newton迭代法求方程
()0
6
cos
2
2x=
-
+
+
=-x
e
x
f x
的根，计算停止的条件为：
6
1
10-
+
<
-
k
k
x
x
；
编程如下：
function newton(f,df,x,a,a0)
syms x
f=input('please enter your equation:') a0=input('please enter you x(0):');
df=diff(f)
e=1e-6;
a1=a0+1;
N=0;
while abs(a1-a0)>e
a=a0-subs(f,a0)/subs(df,a0); a1=a0; a0=a; N=N+1; end
fprintf('a=%0.6f',a) N
运行结果： >> newton
please enter your equation:exp(x)+2^(-x)+2*cos(x)-6 f =exp(x)+2^(-x)+2*cos(x)-6 please enter you x(0):2
df =exp(x)-2^(-x)*log(2)-2*sin(x) a=1.829384 N =4
2．利用Newton 迭代法求多项式
07951.2954.856.104.5x 234=+-+-x x x
的所有实零点，注意重根的问题。

由于不知道重根的个数，采用试探法求重根。

function [X]=newton2() clc;clear; syms x;
delta=1e-06;
f=inline('x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951'); df=inline('4*x^3-16.2*x^2+21.12*x-8.954');
u=inline('(x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951)/(4*x^3-16.2*x^2+21.12*x-8.954)');
du=inline('1-(x^4-27/5*x^3+264/25*x^2-4477/500*x+27951/10000)/(4*x^3-81/5*x^2+528/25*x-4477/500)^2*(12*x^2-162/5*x+528/25)'); j=1;
for i=0:1:3 x0=i;
while 1
x1=x0-feval(u,x0)/feval(du,x0); err=abs(x1-x0); x0=x1;
y=feval(f,x0); if (err<delta) X(j)=x1; j=j+1;
break
end
end
end
X
% 运行结果：
% X =1.1000 1.1000 2.1000 1.1000 四、数值积分
分别用复化梯形公式和Romberg公式计算积分
dx
x
⎰821
的近似值，要求误差不超
过5
10-，并给出节点个数。

1）复化梯形公式
程序代码：
function [I,step]= Trapezoid (f,a,b,eps)
f='1/x';a=2;b=8;eps=1.0e-5;
n=1;
h=(b-a)/2;
I1=2;
I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps
n=n+1;
h=(b-a)/n;
I1=I2;
I2=0;
for i=0:n-1
x=a+h*i;
x1=x+h;
I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+...
subs(sym(f),findsym(sym(f)),x1));
end
end
I=I2
step=n
结果输出
I =
1.38654458753674
step =
53
2）Romberg公式
程序代码：
function approx = romberg (f,a,b,TOL)
clc;clear;
f=inline('1/x');
a=2;b=8;
TOL=1e-05;
A=zeros(20,20);
A(1,1)=(b-a)*(feval(f,a)+feval(f,b))/2;
h=(b-a)/2;
A(2,1)= A(1,1)/2+h*feval(f,a+h);
A(2,2)=(4*A(2,1)-A(1,1))/3;
errest =abs(A(2,2)-A(1,1)/2);
i=2;
while(errest>TOL)
i=i+1;
h =h/2;
sum=0.0;
for j=1:2:2^(i-1)-1
sum=sum+feval(f,a+j*h);
end;
A(i,1)=A(i-1,1)/2+h*sum;
for j=2:i
power=4^(j-1);
A(i,j)=(power*A(i,j-1)-A(i-1,j-1))/(power-1);
end;
errest=abs(A(i,i)-A(i-1,i-1))/2^(i-1);
end;
if (nargout ==0)
s=sprintf('\t\t approximate value of integral: \t %.12f \n',A(i,i));
s=sprintf('%s \t\t error estimate: \t\t\t\t\t %.4e \n',s,errest);
s=sprintf('%s \t\t number of function evaluations: \t %d \n',s,2^(i-1)+1); disp(s)
else
approx=A(i,i);
end
运行结果：
approximate value of integral: 1.386297441871
error estimate: 8.7527e-006
number of function evaluations: 17
五、插值与逼近
1．给定[]1,1-
上的函数
()
2
25
1
1
x
x
f
+
=
，请做如下的插值逼近：
错误！未找到引用源。

构造等距节点分别取5=n ，8=n ，10=n 的Lagrange 插值多项式；
错误！未找到引用源。

构造分段线性取10=n 的Lagrange 插值多项式； 错误！未找到引用源。

取Chebyshev 多项式()()x n x T n arccos cos ⋅=的零点：
πn k x k 21
2cos
-=，n k ,,1 =
作插值节点构造10=n 的插值多项式
()x f 和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。

程序代码
function Lagrange clc;clear;close all; for i=1:3; if i==1 N=5; elseif i==2 N=8; else
N=10; end
f=inline('1/(1+25*x^2)'); x1=zeros(1,N+1); a=-1; b=1;
for i=1:N+1
x1(i)=a+(i-1)*(b-a)/N; y1(i)=feval(f,x1(i)); end syms x ff=0;
for i=1:N+1 f=1;
for j=1:i-1
f=f*(x-x1(j))/(x1(i)-x1(j)); end
for j=i+1:N+1
f=f*(x-x1(j))/(x1(i)-x1(j)); end
ff=f*y1(i)+ff;
f=1;
end
ff=collect(ff,x);
ff=vpa(ff,4);
y=ff;
p=ezplot(y,[a,b]); grid
YLIM([-0.1 0.6]);
if N==5
set(p,'Color','black');
set(p,'LineStyle','--');
lagrange_5=y
elseif N==8
set(p,'Color','r');
set(p,'LineStyle','--');
lagrange_8=y
else
set(p,'Color','b');
set(p,'LineStyle','--')
lagrange_10=y
end
hold on;
xlabel('x');ylabel('y');
title('y=p(x)');hold on;
end
Lag_Cheb();
x=-1:0.01:1;
y=1./(1+25*x.^2);
acu=plot(x,y);grid on;hold on
set(acu,'Color','m');
set(acu,'LineStyle','-');
legend('N=5','N=8','N=10','Cheb,N=10','¾«È·Öµ'); function Lag_Cheb()
f=inline('1/(1+25*x^2)');
N=10;
x1=zeros(1,N+1);
a=-1;
b=1;
for i=1:N+1
x1(i)=cos((2*i-1)*pi/(2*N));
y1(i)=feval(f,x1(i));
end
syms x
ff=0;
for i=1:N+1
f=1;
for j=1:i-1
f=f*(x-x1(j))/(x1(i)-x1(j));
end
for j=i+1:N+1
f=f*(x-x1(j))/(x1(i)-x1(j));
end
ff=f*y1(i)+ff;
f=1;
end
ff=collect(ff,x);
ff=vpa(ff,4);
yy=ff;
ff=collect(ff,x);
yy=ff;
lagrange_chebshev_10=yy
cheb=ezplot(yy,[a,b]); grid on
YLIM([-0.1 0.6]);
set(cheb,'Color','g');
set(cheb,'LineStyle',':');
结果输出
lagrange_5 =.5673+1.202*x^4-1.731*x^2
lagrange_8 =1.+53.69*x^8-102.8*x^6+61.37*x^4-13.20*x^2
lagrange_10 =1.-220.9*x^10+494.9*x^8-381.4*x^6+123.4*x^4-16.86*x^2 lagrange_chebshev_10 =.7413-5.359*x^10-.4e-2*x^9+18.96*x^8-.1321*x^7-
25.78*x^6+.10*x^5+16.81*x^4+.5e-2*x^3-5.336*x^2+.5288e-3*x
2．已知函数值
k x 0 1 2 3 4 5 6 7 8 9 10 k y
2.51
3.30
4.04
4.70
5.22
5.54
5.78
5.40
5.57
5.70
5.80
和边界条件：(0)0.8,(10)0.2s s ''==． 求三次样条插值函数()y s x =并画出其图形．
程序代码：
function Spline3_1(x,y,df0,dfn) format short
x=[0 1 2 3 4 5 6 7 8 9 10];
y=[2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80]; plot(x,y,'g*','MarkerSize',15);hold on; df0=1; dfn=0;
n=length(x); h=zeros(n-1,1); lan=zeros(n-2,1); mu=zeros(n-2,1); g=zeros(n-2,1); m=zeros(n,1);
m(1)=df0;
m(n)=dfn;
for i=1:n-1
h(i)=x(i+1)-x(i);
end
for i=1:n-2
mu(i)=h(i)/(h(i+1)+h(i));
lan(i)=h(i+1)/(h(i+1)+h(i));
g(i)=3*(mu(i)*(y(i+2)-y(i+1))/h(i+1)+lan(i)*(y(i+1)-y(i))/h(i)); end
A=zeros(n-2,n-2);
A(1,1)=2;
A(1,2)=mu(1);
A(n-2,n-2)=lan(n-2);
for i=2:n-2
A(i,i)=2;
A(i,i-1)=mu(i);
A(i-1,i)=lan(i);
end
g(1)=g(1)-lan(1)*df0;
g(n-2)=g(n-2)-mu(n-2)*dfn;
b=A\g;
for i=2:n-1
m(i)=b(i-1);
end
syms z
for i=1:n-1
xx=x(i):0.01:x(i+1);
sx1=y(i)*(xx-x(i+1)).^2.*(h(i)+2*(xx-x(i)))/h(i)^3;
sx2=y(i+1)*(xx-x(i)).^2.*(h(i)-2*(xx-x(i+1)))/h(i)^3;
sx3=m(i)*(xx-x(i+1)).^2.*(xx-x(i))/h(i)^2;
sx4=m(i+1)*(xx-x(i)).^2.*(xx-x(i+1))/h(i)^2;
sx=sx1+sx2+sx3+sx4;
z1=y(i)*(z-x(i+1)).^2.*(h(i)+2*(z-x(i)))/h(i)^3;
z2=y(i+1)*(z-x(i)).^2.*(h(i)-2*(z-x(i+1)))/h(i)^3;
z3=m(i)*(z-x(i+1)).^2.*(z-x(i))/h(i)^2;
z4=m(i+1)*(z-x(i)).^2.*(z-x(i+1))/h(i)^2;
zz=z1+z2+z3+z4;
zz= collect(zz);
vpa(zz,4)
plot(xx,sx,'r');
hold on;
grid on;
end
程序输出
ans =2.510+.1379*z^3-.3479*z^2+z
ans =2.692-.4369e-1*z^3+.1969*z^2+.4552*z ans =2.287+.6872e-2*z^3-.1065*z^2+1.062*z ans =3.655-.4380e-1*z^3+.3495*z^2-.3060*z ans =-6.080+.1083*z^3-1.476*z^2+6.995*z ans =41.14-.2694*z^3+4.190*z^2-21.34*z
ans =-109.8+.4295*z^3-8.390*z^2+54.14*z ans =133.0-.2784*z^3+6.475*z^2-49.91*z
ans =-57.75+.9411e-1*z^3-2.465*z^2+21.61*z ans =75.03-.8804e-1*z^3+2.453*z^2-22.65*z。