vpa function did not convert sym results as intended

1 view (last 30 days)
I have a function that applied gaussian elimination to a matrix and is suppose to eliminate values from a few rows and columns. A few days ago, the code was working fine and I got a bunch of complex elements as the output of the user defined function below. Today, the output is just a sym format and not double complex. I cannot figure out what the problem is.
I have attached the csv file of the ybus matrix.
function yreduced=gauss_elimination(ybus,yb,Nb)
indexyb=ybus(1,:);
MM=ybus(2:end,2:end); % Augmented Matrix
MM=sym(MM);
%%%%%%%%%%%%%%%% Gauss elimination method %%%%%%%%%%%%%%
disp('Gauss elimination method:');
[m]=size(MM);
y1=size(Nb);
y2=sort(Nb,'descend')
for j=1:y1
for k=1:m
if k==y2(j)
for z=k-3:m
if MM(k,k)==0
disp("replace row"+k + "with row" +z)
tmp=MM(k,:); MM(k,:)=MM(z,:);
MM(z,:)=tmp;
end
end
% end
i=k-3;
ss=1;
while ss<=3
disp("Eliminate ybus("+i+","+k+")")
MM(i,:)=MM(i,:)-MM(k,:)*(MM(i,k)/MM(k,k))
% i=i-3;
ss=ss+3;
end
end
end
end
% yreduced=M;
yreduced=vpa(MM,6);
  2 Comments
Torsten
Torsten on 30 May 2022
Why do you define MM as sym ?
MM is a numerical matrix from the .csv file. You will not gain higher precision this way.
Leonie Bule
Leonie Bule on 30 May 2022
Edited: Leonie Bule on 30 May 2022
I have a live script so I declared MM as sym for algebraic visualization. I have noticed a problem with the for loop
for z=k-3:m
if MM(k,k)==0
disp("replace row"+k + "with row" +z)
tmp=MM(k,:); MM(k,:)=MM(z,:);
MM(z,:)=tmp;
end
end
It doesn't replace the MM(k,k) k rows with non-zero rows. Maybe the sym and vpa function fail when the elements are not double precision?

Sign in to comment.

Answers (1)

SAI SRUJAN
SAI SRUJAN on 27 Oct 2023
Hi Leonie Bule,
I understand that you are facing an issue regarding the implementation of gauss elimination method.
You can follow the given example to proceed further.
a = [3 4 -2 2
4 9 -3 5
-2 -3 7 6
1 4 6 7 ];
a=sym(a);
[m,n]=size(a);
for j=1:m-1
for z=2:m
if a(j,j)==0
t=a(j,:);
a(j,:)=a(z,:);
a(z,:)=t;
end
end
for i=j+1:m
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
end
end
vpa(a,6);
You can refer to the following documentation to understand more about "vpa" MATLAB function.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!