Gaussian Elimination in Matlab - Programmers Heaven

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Categories

Gaussian Elimination in Matlab

veronicak5678veronicak5678 Posts: 4Member
I am trying to write a Matlab program to perform Gaussian Elimination with Partial Pivoting. I getting bad results, but I can't see where my code deviates from the algorithm I'm using. Here is a link to the book with the steps I'm trying to follow:

[link=http://books.google.com/books?id=wmcL0y2avuUC&pg=PA368&lpg=PA368&dq=gaussian+elimination+algorithm+6.2+numerical+analysis&source=bl&ots=1esN2ZdIPZ&sig=jQkt8_9lM8CloVzy3vWz8vVxoR4&hl=en&ei=RvssTd_UD42usAOZ9ImbBw&sa=X&oi=book_result&ct=result&resnum=1&ved=0CBcQ6AEwAA#v=onepage&q&f=false]http://books.google.com/books?id=wmcL0y2avuUC&pg=PA368&lpg=PA368&dq=gaussian+elimination+algorithm+6.2+numerical+analysis&source=bl&ots=1esN2ZdIPZ&sig=jQkt8_9lM8CloVzy3vWz8vVxoR4&hl=en&ei=RvssTd_UD42usAOZ9ImbBw&sa=X&oi=book_result&ct=result&resnum=1&ved=0CBcQ6AEwAA#v=onepage&q&f=false[/link]

This is algorithm 6.2 which begins on page 362. Here is my code:
[code]clear
A= [1 -1 2 -1 ; 2 -2 3 -3; 1 1 1 0; 1 -1 4 3 ]
[~,q]=size(A);
n=q-1;
x=zeros(n,1);
m=[n,q];
p=0;
nrow=ones(n,1);
ncopy = [1, q];
%-------------------------------Step 1
for i=1:n
nrow(i,1)=i;
end
%------------------------------Step 2
for i=1:n-1
%-----------------------------Step 3
maxval = 0;
index=0;
for j =i:n
if abs(A(nrow(j),i))>maxval
maxval=abs(A(nrow(j),i));
index = j;
end
end
p=index;
%-----------------------------Step 4
if A(nrow(p),i) ==0
display 'No!'
end
%-----------------------------Step 5
if nrow(i,:)~=nrow(p,:)
ncopy = nrow(i,:);
nrow(i,:)=nrow(p,:);
nrow(p,:)= ncopy;
end
%------------------------------Step 6
for j=i+1:n
%------------------------------Step 7
m(nrow(j),i) =A(nrow(j),i)/A(nrow(i),i);
%------------------------------Step 8
A(nrow(j),:) =A(nrow(j),:)- (m(nrow(j),i)*A(nrow(i),:));
end
end
%------------------------------Step 9
if A(nrow(n),n)==0
display 'No!'
end
%------------------------------Step 10
x(n,1)=A(nrow(n),q)/A(nrow(n),n);
%------------------------------Step 11
for i = n-1:-1:1
sig=0;
for j=i+1:n
sig =sig + A(nrow(i),j)*x(j,1);
end
x(i,1) = (A(nrow(i),q) - sig)/A(nrow(i),i);
end
%-----------------------------Step 10
display ' A is'
disp (A)
display ' x is'
disp (x)[/code]
Sign In or Register to comment.