# Finding a particular solution in Matlab

So I have a question about finding the particular solution to a non-homogeneous system of equations in Matlab. This is building on an assignment where I had to find the spanning set (I think this is the same as the column space?) of a system. The snippet of code shows my solution to the problem and here is my reasoning:

Begin by finding the ordered set of indexes of the pivot and the non-pivot positions in the reduced coefficient matrix A (with dimensions m x n).

Next, reduce the augmented matrix and call it R.

The first for loop sets all of the free variables (non-pivot positions) in R to zero. All I have left are the basic variables and the solution column. Since all pivot positions are equal to 1, I can say that variable is equal to the solution column plus some combination of free variables. (Right now, I am only interested in the solution column - for example, if x1 = 3 + x2, I only want the 3)

The second for loop assigns entries in the last column to a row in matrix p for the particular solution.

My code has worked for all the test cases I have tried so far, but I am wondering whether my solution is actually correct or if I just got lucky on a few examples? If I'm wrong (or if I'm right) is there a better way to solve this?

%finds the ordered set of indexes of the pivot columns
[~,pivot_c]=rref(A);

%finds the ordered set of indexes of the non-pivot columns
S=1:n;
nonpivot_c=setdiff(S,pivot_c);

%set R equal to the reduced echelon form of the augmented matrix
R = rref([A,b]);

%pre-allocate space for the particular solution p
p = zeros(n,1);

%set all free variables in R equal to zero
for i = 1:n
if i == nonpivot_c
R(:,i) = zeros(m,1);
end
end

%now the only values remaining are the basic variables
%the last column is the values of the particular solution
%when the free variables are equal to zero
for i = 1:m
for j = 1:n
if R(i,j) == 1
p(j,1) = R(i, (n+1));
end
end
end

%display the solution
p


• I would suggest computational science SE or Stack Overflow. – K.defaoite Feb 23 at 3:12

What you have done usually works (as long as it is known a priori that the equation $$Ax = b$$ has a solution). Unfortunately, if i == nonpivot_c does not do what you think it does: Matlab will never treat i == nonpivot_c as a true statement (and will therefore never zero out the free columns) because it is an array of Booleans. If any of your free columns happens to contain $$1$$ as an entry, this will lead to trouble in your final for loop.

To correctly check whether i is an element of nonpivot_c, use if ismember(i,nonpivot_c).

In any case, the method you outlined is not the most efficient way to do what you're trying to do. Consider the following approach:

A = [0,0,2,3,4;...
1,2,3,4,5;...
-1,-2,-5,-7,-9];
b = [2;0;-2];        % Has a solution
% b = [1;0;0];         % Has no solution

[m,n] = size(A);
[R,p_cols] = rref([A,b]);
if p_cols(end)==n+1
disp('---NO SOLUTION---')
else
p=zeros(n,1);
for i = 1:length(p_cols)
p(p_cols(i)) = R(i,end);
end
disp('SOLUTION:')
disp(p)
end


As a bonus, we can eliminate the for loop:

[m,n] = size(A);
[R,p_cols] = rref([A,b]);
if p_cols(end)==n+1
disp('---NO SOLUTION---')
else
p=zeros(n,1);
p(p_cols) = R(1:length(p_cols),end);
disp('SOLUTION:')
disp(p)
end

• Thanks, that helped a lot. Also way more concise than what I had. – StupidQuestionsGuy Feb 23 at 14:18