I'm currently trying to implement a finite element method in Python and I have to solve a linear system where the matrix is under determined, but the system should be well-defined by the boundary condition.
To be precise, in a very simple case, the matrix looks like this
[ 12., 6., -12., 6.],
[ 6., 4., -6., 2.],
[-12., -6., 12., -6.],
[ 6., 2., -6., 4.]]
And the right-hand side is just an array of numbers.
Obviously, the third column is the negative of the first, so the matrix is under determined.
However, I want to impose that the first entry in the solution vector is zero, so the solution should be of the form [0,x,y,z]. Now we have 3 unknowns and 3 equations, so there should be a unique solution, but I haven't found a scipy function that does that, since the normal scipy.linalg.solve only takes the matrix and the right-hand side as inputs.
I've tried to rewrite the linear system, but in some cases there are undefined variables in the right-hand side as well and more set values in the unknown vector, so I'm looking for a procedure that solves an any linear system with the matrix, and two vectors with unknowns.
You're question is a bit confusing because your matrix has rank 2 and not 3 (this means it gives you only 2 linearly independent equations). Now if I call your matrix A and the right hand side b the best you can hope for is a solution of the form
A@[0,0,x,y]=b.
Keep in mind that this does not always work. Indeed since the last two columns of A are linearly independent this solution exists if and only if b is in the image of A. So b must come from some multiplication of A with a vector.
Now if we remember one simple fact from linear algebra we can find you a solution in the form you hope for. We can write any solution x to an underdetermined linear equation as x=v+n with v being any solution and n being from the null space of the matrix. So we just find a solution that is not in the form yet and then bring it into the form by moving it threw the null space of A.
Let's go thew an explicit example: suppose we have
b=[12,4,-12,8](=A@[1,2,3,4]).
(Of course I assume here we don't know that b=A@[1,2,3,4].)
We can find a solution with v=np.linalg.lstsq(A, b)[0](=[ 0.4, -0.8, -0.4, 1.2]). Here we should double check that A@v==b (or at least that np.isclose(A@v,b)) since np.linalg.lstsq does not throw an exception if no solution exists but gives you an approximate one.
If now you want a different solution in your special form you just calculate the null space
null_space = scipy.linalg.null_space(A)
and then how we have to change v to get it to be zero in the first two entries
coefs = np.linalg.solve(null_space[:2], v[:2])
w = -null_space@coefs
x = v+w.
Turns out x=[0,0,0,2].