以下代码是否可以避免这些舍入错误?
Are these round-off errors avoidable in the following code?
我正在尝试在 Matlab 中编写一个程序来使用 LU 分解求解线性方程组,该分解使用高斯消元法,因此需要很多算术步骤。
答案接近于正确的解决方案,但与 Python 等其他语言相比,舍入误差相当高
例如,其中一个解正好是 3,但我得到 2.9877。
我知道内置函数应该用于这些琐碎的事情,因为 Matlab 是一种高性能计算语言,但如果我仍然想用循环等来做,我总是会遇到舍入错误或者有什么办法在进行数值计算时减少这些?
我附上了代码,但它很大,不值得一读。为了完成,我仍然附上它。
可以注意到许多算术运算的使用会引入很多舍入误差。
这些舍入误差是 Matlab 固有的并且不可避免吗?
clc
clear
%No of equations is n
n=3;
%WRITING THE Coefficients
A(1,1)=3;
A(1,2)=-0.1;
A(1,3)=-0.2;
B(1)=7.85;
A(2,1)=0.1;
A(2,2)=7;
A(2,3)=-0.3;
B(2)=-19.3;
A(3,1)=0.3;
A(3,2)=-0.2;
A(3,3)=10;
B(3)=71.4;
%Forward Elimination
for i=1:n-1
for j=i+1:n
fact=A(j,i)/A(i,i);
A(j,i)=fact;
A(j,j:n)=A(j,j:n)-fact*A(i,j:n);
B(j)=B(j)-fact*B(i);
end
end
disp(A)
% Calculating d matrix
sum=0;
D(1)=B(1);
for i=2:n
for j=1:i-1
sum=sum+A(i,j)*B(j);
D(i)=B(i)-sum;
end
end
disp("D =")
disp(transpose(D))
%Back Substitution
X(n)=D(n)/A(n,n);
for z=n-1:-1:1
sum=0;
for w=z+1:n
sum=sum+A(z,w)*X(w);
end
X(z)=(D(z)-sum)/(A(z,z));
end
disp("Solution X is:")
disp(transpose(X))
永远不要忘记不信任你的编码。
如果您注释掉第 24 行并将第 26 行替换为 A(j,i:n)=A(j,i:n)-fact*A(i,i:n);
,您将得到很好的解决方案(长格式)
3.000000000000000
-2.500000000000000
7.000000000000002
我并不是说这是最好的修复(它不是),但它清楚地表明舍入是无罪的。错误的解决方案仍然很接近,因为系统强烈 diagonal-dominant.
我正在尝试在 Matlab 中编写一个程序来使用 LU 分解求解线性方程组,该分解使用高斯消元法,因此需要很多算术步骤。
答案接近于正确的解决方案,但与 Python 等其他语言相比,舍入误差相当高
例如,其中一个解正好是 3,但我得到 2.9877。
我知道内置函数应该用于这些琐碎的事情,因为 Matlab 是一种高性能计算语言,但如果我仍然想用循环等来做,我总是会遇到舍入错误或者有什么办法在进行数值计算时减少这些?
我附上了代码,但它很大,不值得一读。为了完成,我仍然附上它。 可以注意到许多算术运算的使用会引入很多舍入误差。
这些舍入误差是 Matlab 固有的并且不可避免吗?
clc
clear
%No of equations is n
n=3;
%WRITING THE Coefficients
A(1,1)=3;
A(1,2)=-0.1;
A(1,3)=-0.2;
B(1)=7.85;
A(2,1)=0.1;
A(2,2)=7;
A(2,3)=-0.3;
B(2)=-19.3;
A(3,1)=0.3;
A(3,2)=-0.2;
A(3,3)=10;
B(3)=71.4;
%Forward Elimination
for i=1:n-1
for j=i+1:n
fact=A(j,i)/A(i,i);
A(j,i)=fact;
A(j,j:n)=A(j,j:n)-fact*A(i,j:n);
B(j)=B(j)-fact*B(i);
end
end
disp(A)
% Calculating d matrix
sum=0;
D(1)=B(1);
for i=2:n
for j=1:i-1
sum=sum+A(i,j)*B(j);
D(i)=B(i)-sum;
end
end
disp("D =")
disp(transpose(D))
%Back Substitution
X(n)=D(n)/A(n,n);
for z=n-1:-1:1
sum=0;
for w=z+1:n
sum=sum+A(z,w)*X(w);
end
X(z)=(D(z)-sum)/(A(z,z));
end
disp("Solution X is:")
disp(transpose(X))
永远不要忘记不信任你的编码。
如果您注释掉第 24 行并将第 26 行替换为 A(j,i:n)=A(j,i:n)-fact*A(i,i:n);
,您将得到很好的解决方案(长格式)
3.000000000000000
-2.500000000000000
7.000000000000002
我并不是说这是最好的修复(它不是),但它清楚地表明舍入是无罪的。错误的解决方案仍然很接近,因为系统强烈 diagonal-dominant.