Is this algorithm correct?
2 views (last 30 days)
Show older comments
%% 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 Comments
Answers (1)
Jan
on 23 Apr 2014
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 .
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!