Is this algorithm correct?
显示 更早的评论
%% Part 1: Crout factorization
n=10 A = full(gallery('tridiag',n,-1,2,-1)) i = 2:(n-1); bmid = i.^2 / ((n+1).^4) b = [1+1/(n+1)^4, bmid, 6+(n^2)/(n+1)^4]' for i = 1:n L(i,1) = A(i,1) end
for j = 1:n U(1,j) = A(1,j)/L(1,1) end
for j = 2:n for i = j:n sum = 0.0 for k = 1:(j-1) sum = sum + L(i,k) * U(k,j) end L(i,j) = A(i,j) - sum end
U(j,j) = 1;
for i = (j+1):n
sum = 0.0
for k = 1:(j-1)
sum = sum + L(j,k) * U(k,i)
end
U(j,i) = (A(j,i) - sum)/L(j,j)
end
end
x=A\b
error=max(abs(b-A*x))
%% Part 2: Gauss Seidel iteration n = 10 n = 10 A = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1); L = 2*diag(ones(n,1)) - diag(ones(n-1,1),-1); U = -diag(ones(n-1,1),1);
b = 1/(n+1)^4*(1:n)'.^2
b(1) = b(1) + 1; b(n) = b(n) + 6
x = zeros(n,1)
X = ones(n,1)
tol = 10e-8;
epsilon = 10*tol
while (epsilon > tol)
x = inv(L)*(b-U*x)
epsilon = max(abs(X-x))
X = x;
end
Can someone tell me are these algorithms correct or not because when I try to perform 1000x1000 matrix, it is taking me more than a hour to compute the result.
回答(1 个)
Jan
2014-4-23
0 个投票
You ask for correctness but argument with the speed.
Obviously your code is slow and there are many ways to improve it substantially. E.g. using Matlab's sum function is much faster than overwrite this important function by a local variable and sum in a FOR loop.
The multiplication with the inverse is a big DON'T in numerical software. See help slash .
类别
在 帮助中心 和 File Exchange 中查找有关 Image Arithmetic 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!