大连理工大学矩阵分析matlab上机作业

x=zeros(n,1); %为列向量 x 预分配存储空间 y=1:n; %定义行向量 y y=y'; %把行向量 y 改成列向量 for i=1:n
x(i)=1/i; %按要求给向量 x 赋值，其值递减 end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 2 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的无穷范数 normxinf normy1=norm(y,1); %求解向量 y 的 1 范数 normy1 normy2=norm(y,2); %求解向量 y 的 2 范数 normy2 normyinf=norm(y,inf); %求解向量 y 的无穷范数 normyinf z1=[normx1,normx2,normxinf]; z2=[normy1,normy2,normyinf]; end
for i=2:n
for j=i:n U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
式
%Doolittle 分解计算上三角矩阵的公
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); %Doolittle 分解计算下三角矩 阵的公式
end
1 1 1 ������ x = (1, 2 , 3 , … , ������) ,
������ = (1,2, … , ������)������.
对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改
你的程序使得计算结果相对精确呢？
1.1 源代码
function [z1,z2]=norm_vector(n) %向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三 种范数，n 为输入参数
矩阵与数值分析
上机实验报告
学生姓名：
高明星
班 级： 矩阵与数值分析 4 班
学 号：
31603026
专 业： 船舶与海洋工程
授课老师：
董波
大连理工大学
Dalian University of Technology
1 计算给定向量的范数
题目：考虑计算给定向量的范数：输入向量x = (������1, ������2, … , ������������)������,输出 ‖������‖1, ‖������‖2, ‖������‖∞。请编制一个通用程序，并用你编制的程序计算如下向量的范 数：
1000.0000000000
结果分析：该源代码在求解向量 x 的 1 范数和 2 范数时存在“大数吃小数”的问
题，会降低其精度，带来误差，当 n 越大时，误差会越大，因此可逆序求和。修
改后的程序代码为:
function z1=norm_vector2(n)
%向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三种范数
y=lnd/(d-1) end if 用此算法计算x ∈ [−10−15, 10−15]时的函数值，画出图像。比较一下发生了什 么？
2.1 源程序及运行结果
clc x=(-1:0.01:1)*1E-15;
y1=log(1+x)./x; d=1+x; for i=1:length(x)
if d(i)==1 y2(i)=1;
n=
1000
z1 =
7.485470860550344 1.282160117411847 1.000000000000000
z2 =
1.0e+05 * 5.005000000000000 0.182711110773264 将以上数据汇总为表 1-1
0.010000000000000
表 1-1 向量 x,y 各范数汇总
3 秦九韶算法
题目：首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用 程序，你的程序包括输入多项式的系数以及给定点，输出函数值。利用你编制 的程序计算
p(x) = (������ − 2)9 = ������9 − 18������8 + 144������7 − 672������6 + 2016������5 − 4032������4 + 5376������3 − 4608������2 − 512
用该程序求得的解与真实解绘制的对比图，可以发现，在零点领域内，用秦九韶
求得的解是不稳定的，这主要是由于相近数相减损失了有效数字 。
4 LU 分解和 PLU 分解
题目：编制计算给定矩阵A的LU分解和PLU分解的通用程序，然后用你编制的程 序完成下面两个计算任务：
(1) 考虑
1
−1 A= ⋮
−1 [−1
%n 为输入参数
x=zeros(n,1); %为列向量 x 预分配存储空间
for i=1:n x(n-i+1)=1/(i); %按升序给向量 x 赋值，其值递增
end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 1 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的 1 范数 normxinf
0⋯0 1 ⋱⁡⁡ ⋱ ⁡⁡⁡⁡⁡⁡ ⋮ ⁡⁡⁡⁡⁡⁡⁡ ⋮ ⋱ ⋱ 0 1 ∈ ������������×������, ⋯ −1⁡⁡⁡⁡⁡⁡⁡1 1 ⋯ −1⁡⁡⁡⁡−1 1]
自己取定x ∈ ������������，并计算b = Ax。然后用你编制的不选主元和列主元的Gauss
消去法求解该方程组，记你计算出的解为���̂���。对n从5到30估计计算解的精
U(i,n)=1;
end
A=L+D+U;
end LU 分解函数 LU (A)定义：
function [L,U] =LU (A) %求矩阵 A 的 LU 分解，返回单位下三角矩阵 L，上三角矩阵 U
n=length(A);
L=eye(n);
U=zeros(n);
U(1,:)=A(1,:);
L(:,1)=A(:,1)/U(1,1);
在x=2邻域附近的值。画出p(x)在x ∈ [1.95,20.5]上的图像。
3.1 源程序
用秦九韶算法计算一个多项式在给定点的函数值的通用程序为： function y= QJS(a,x)
%QJS()为用秦九韶算法计算多项式值的函数 %a 为多项式的系数向量，x 为自变量，y 为因变量 n=length(a); y=a(1); for i=2:n y=y*x+a(i); end end
end
end PLU 分解函数 PLU (A)定义：
function [P,L,U] =PLU (A) %求矩阵 A 的 LU 分解，返回单位下三角矩阵 L，上三角矩阵 U，置换矩阵 P
n=length(A);
P=eye(n);
L=eye(n); index=1; %用于显示列主元所在行号
for j=1:n-1
616
599
7
0
2 计算函数的值
题目：考虑y = f(x) = ln⁡(1+������),其中定义f(0)=1，此时f(x)是连续函数。用此公
������
式计算当x ∈ [−10−15, 10−15]时的函数值，画出图像。另一方面，考虑下面算 法：
d=1+x if d=1 then
y=1 else
结果分析：图 2-1 为采用两种算法计算函数 y=ln(1+x)/x 在区间[10-15,1015]上的
图像对比，可以看出第一种算法所画图像在 y=1 处上下震荡，呈现出数值的不稳
定，这主要是由于小数做除数损失了有效数字，而第二种算法避免了有效数字的
损失，所画出的图像在 x=0 附近均为 1，数值很稳定。
n=10
n=100
n=1000
x
1
2.928968253968254
5.187377517639621
7.485470860550344
x
2
1.244896674895769
1.278664889713052
1.282160117411847
x

1.000000000000000
1.000000000000000
y= 0
1.5 x 10-11 1
秦九韶公式 y=(x-2)9
0.5
y
0
-0.5
-1
-1.5
1.95
2
2.05
x
图 3-1 y=(x-2)9 的图像
结果分析：当 x=2 时，通过调用秦九韶函数 QJS()计算多项式 y=(x-2)9，得到 y=0，
初步验证了该通用程序的正确性。图 3-1 为多项式 y=(x-2)9 在区间[1.95,20.5]上
55.000000000000000 19.621416870348583 10.000000000000000 n=
100
z1 =
5.187377517639621 1.278664889713052 1.000000000000000 z2 =
1.0e+03 *
5.050000000000000 0.581678605417115 0.100000000000000
else y2(i)=log(d(i))/(d(i)-1);
end end plot(x,y1,'--r',x,y2,'m'); xlabel('x'); ylabel('y'); legend('修改前','修改后');