Finding a particular solution in Matlab

Question

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
 

Answer 1

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