Matlab debugging for LU decomposition

  • Thread starter joehardy
  • Start date
  • #1
1
0
these are my mfiles for Lu decomposition with partial pivoting and substitution

1.LUdecomp.m
function [x]=LUdecomp(a,b,n,tol,er,x)
for i=1:n
o(i)=0;
s(i)=0;
end
er=0;
decompose(a,n,tol,o,s,er)
if er~=-1
substitute(a,o,n,b,x);
end
end

2.decompose.m
function [a]=decompose(a,n,tol,o,s,er)
for i=1:n
o(i)=i;
s(i)=abs(a(i,1));
for j=2:n
if abs(a(i,j))>s(i)
s(i)=abs(a(i,j));
end
end
end
for k=1:n-1
pivot(a,o,s,n,k)
if abs(a(o(k),k)/s(o(k)))<tol
er=-1;
abs(a(o(k),k)/s(o(k)));
break
end
for i=k+1:n
factor=a(o(i),k)/a(o(k),k);
a(o(i),k)=factor;
for j=k+1:n
a(o(i),j)=a(o(i),j)-factor*a(o(k),j);
end
end
end
if abs(a(o(k),k)/s(o(k)))<tol
er=-1;
abs(a(o(k),k)/s(o(k)));
end
end

3.pivot.m
function [o]=pivot(a,o,s,n,k)
p=k;
big=abs(a(o(k),k)/s(o(k)));
for i=k+1:n
dummy=abs(a(o(i),k)/s(o(i)));
if dummy>big
big=dummy;
p=i;
end
end
dummy=o(p);
o(p)=o(k);
o(k)=dummy;
end

4.substitute.m

function [x]=substitute(a,o,n,b,x)
for i=2:n
sum=b(o(i));
for j=1:i-1
sum=sum-a(o(i),j)*b(o(j));
end
b(o(i))=sum;
end
x(n)=b(o(n))/a(o(n),n);
for i=n-1:-1:1
sum=0;
for j=i+1:n
sum=sum+a(o(i),j)*x(j);
end
x(i)=(b(o(i))-sum)/a(o(i),i);
end
end
thank you
 

Answers and Replies

  • #2
34,866
6,598
What's your question? From the thread title, you seem to be looking for help in debugging this code. Can you narrow down which function is not working correctly?

Also, it's a good idea to surround your code with [ code] and [ /code] tags (omit the leading space.
these are my mfiles for Lu decomposition with partial pivoting and substitution

1.LUdecomp.m
Code:
function [x]=LUdecomp(a,b,n,tol,er,x)
for i=1:n
    o(i)=0;
    s(i)=0;
end
er=0;
decompose(a,n,tol,o,s,er)
if er~=-1
    substitute(a,o,n,b,x);
end
end
2.decompose.m
Code:
function [a]=decompose(a,n,tol,o,s,er)
    for i=1:n
        o(i)=i;
        s(i)=abs(a(i,1));
        for j=2:n
            if abs(a(i,j))>s(i)
            s(i)=abs(a(i,j)); 
            end
        end
    end
for k=1:n-1
    pivot(a,o,s,n,k)
    if abs(a(o(k),k)/s(o(k)))<tol
        er=-1;
        abs(a(o(k),k)/s(o(k)));
    break
    end
    for i=k+1:n
        factor=a(o(i),k)/a(o(k),k);
        a(o(i),k)=factor;
        for j=k+1:n
            a(o(i),j)=a(o(i),j)-factor*a(o(k),j);
        end
    end
end
if abs(a(o(k),k)/s(o(k)))<tol
        er=-1;
        abs(a(o(k),k)/s(o(k)));
end
end

3.pivot.m
Code:
function [o]=pivot(a,o,s,n,k)
p=k;
big=abs(a(o(k),k)/s(o(k)));
for i=k+1:n
    dummy=abs(a(o(i),k)/s(o(i)));
    if dummy>big
        big=dummy;
        p=i;
    end
end
dummy=o(p);
o(p)=o(k);
o(k)=dummy;
end
4.substitute.m
Code:
function [x]=substitute(a,o,n,b,x)
for i=2:n
    sum=b(o(i));
    for j=1:i-1
        sum=sum-a(o(i),j)*b(o(j));
    end
    b(o(i))=sum;
end
x(n)=b(o(n))/a(o(n),n);
for i=n-1:-1:1
    sum=0;
    for j=i+1:n
        sum=sum+a(o(i),j)*x(j);
    end
    x(i)=(b(o(i))-sum)/a(o(i),i);
end
end
thank you
 

Related Threads on Matlab debugging for LU decomposition

  • Last Post
Replies
6
Views
5K
  • Last Post
Replies
1
Views
966
Replies
7
Views
6K
  • Last Post
Replies
3
Views
1K
Replies
1
Views
2K
  • Last Post
Replies
6
Views
1K
  • Last Post
Replies
6
Views
5K
Replies
2
Views
900
Replies
16
Views
8K
  • Last Post
Replies
0
Views
1K
Top