Is this algorithm correct?
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
%% 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 Commenti
Risposte (1)
Jan
il 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 .
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!