计算方法上机题.docx

计算,要求误差小于,给出实现算法。
算法组织:利用绝对值的递增相加可以减小误差
具体算法:
%% 公式求和%%
S = 0;
k = 100000;
%%利用for循环对公式进行求和
for n = 1: 100000
S = S+1/(k*k);
k = k-1;
end
%%精确到小数点后6位,以保证误差范围
vpa(S,7)
运行程序后,,精确到了小数点后6位。
编写实现对N阶非奇矩阵A进行LU分解的程序。
算法组织:
根据定理可知若N阶矩阵的各阶顺序主子式行列式不为零,则存在唯一的单位下三角矩阵L和上三角矩阵U,满足A=LU,因此首先需要判断A中的元素akk(k=0,1,2…,n)不为零。
根据如下的式子计算L和U中的元素
当i=1时,得到u1j=a1j(j=1,2,3,…,n),而当j=1时,得li1= ai1/u11,由此可以计算出L的第一列元素和U的第一行元素。然后可以依据下式计算L和U中的其他元素。
具体算法:
[n,n]=size(A); %检测A的阶数
L=eye(n); %初始化单位下三角阵L
U=zeros(n,n); %初始化上三角阵U
for i=1:n
U(1,i)=A(1,i); %给上三角阵U的第一行赋值
L(i,1)=A(i,1)/U(1,1); %给单位下三角阵L的第一列赋值
end
for i=2:n
for j=i:n
for k=1:i-1
M(k)=L(i,k)*U(k,j);
end
U(i,j)=A(i,j)-sum(M); %给上三角阵U的第i行赋值
end
for j=i+1:n
for k=1:i-1
M(k)=L(j,k)*U(k,i);
end
L(j,i)=(A(j,i)-sum(M))/U(i,i); %给单位下三角阵L的第i列赋值
end
end
%分别输出L,U的矩阵
L
U
,并对所附的方程组进行求解。具体的要求参见所附的相关文档。
针对本专业中所碰到的实际问题,提炼一个使用方程组进行求解的例子,并对求解过程进行分析、求解。
具体算法:
fid = fopen('','r'); % 获取数据文件文件头
headfile = fread(fid,3,'uint32'); % 确认数据文件标识号是否正确
id = headfile(1);
if id ~= hex2dec('F1E1D1A0')
fprintf('不是正确的数据文件.\n');
return;
end
% 获取数据文件版本号
ver = headfile(2);
if ver == hex2dec('101')
% 获取方程阶数
n = headfile(3)
% 获取稀疏矩阵
A = fread(fid,[n,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n);
end
if ver == hex2dec('102')
% 系数矩阵为条状对角阵
% 获取方程阶数
n = headfile(3)
% 获取稀疏矩阵
A = fread(fid,[n,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n);
end
if ver == hex2dec('201')
% 系数矩阵为压缩格式下的条状对角阵
% 获取方程阶数
n = headfile(3)
% 获取条状矩阵的上下带宽
headfile2 = fread(fid,2,'uint32');
p = headfile2(1);
q = headfile2(2);
% 压缩系数矩阵
A = fread(fid,[p+q+1,n],'float32');
% 获取右端系数
B = fread(fid,[n,1],'float32');
% 对系数矩阵A和右端向量的方程组使用高斯方法求解
x=GAUSSPP(A,B,n,p,q);
end
fclose(fid);
GAUSSPP(A,B,n)的程序如下:
function x=GAUSSPP(A,B,n)
for k=1:n-1
for i=k+1:n
l(