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.

Answer 1

[Jan]

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 .

[Sam]

Yes because my instructor said if my algorithm is wrong, then it can take long time.