Is this algorithm correct?
2 次查看(过去 30 天)
显示 更早的评论
%% 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.
0 个评论
回答(1 个)
Jan
2014-4-23
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 .
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!