• Support PF! Buy your school textbooks, materials and every day products Here!

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
33,508
5,193
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
910
Replies
7
Views
6K
  • Last Post
Replies
3
Views
1K
Replies
1
Views
1K
  • Last Post
Replies
6
Views
1K
  • Last Post
Replies
6
Views
4K
Replies
16
Views
7K
Replies
2
Views
847
  • Last Post
Replies
0
Views
1K
Top