solve() only works for nonsingular systems of equations.

Use a generalized inverse for singular systems:

> A<- matrix(c(1,2,1,1, 3,0,0,4, 1,-4,-2,-2, 0,0,0,0), ncol=4, byrow=TRUE)

> A

[,1] [,2] [,3] [,4]

[1,] 1 2 1 1

[2,] 3 0 0 4

[3,] 1 -4 -2 -2

[4,] 0 0 0 0

> b<- c(0,2,2,0) #rhs

> b

[1] 0 2 2 0

>

> require('MASS')

> giA<- ginv(A) #M-P generalized inverse

> giA

[,1] [,2] [,3] [,4]

[1,] 0.6666667 1.431553e-16 0.33333333 0

[2,] 0.3333333 -1.000000e-01 -0.03333333 0

[3,] 0.1666667 -5.000000e-02 -0.01666667 0

[4,] -0.5000000 2.500000e-01 -0.25000000 0

>

> require('Matrix')

> I<- as.matrix(Diagonal(4)) #order 4 identity matrix

> I

[,1] [,2] [,3] [,4]

[1,] 1 0 0 0

[2,] 0 1 0 0

[3,] 0 0 1 0

[4,] 0 0 0 1

>

> giA%*%b #particular solution

[,1]

[1,] 6.666667e-01

[2,] -2.666667e-01

[3,] -1.333333e-01

[4,] -2.220446e-16

> giA%*%A - I #matrix for parametric homogeneous solution

[,1] [,2] [,3] [,4]

[1,] 0.000000e+00 0.0 0.0 5.551115e-16

[2,] 3.469447e-17 -0.2 0.4 4.024558e-16

[3,] 4.510281e-17 0.4 -0.8 2.706169e-16

[4,] -3.330669e-16 0.0 0.0 -7.771561e-16

At 09:34 PM 5/21/2011, dslowik wrote:

>I have a simple system of linear equations to solve for X, aX=b:

> > a

> [,1] [,2] [,3] [,4]

>[1,] 1 2 1 1

>[2,] 3 0 0 4

>[3,] 1 -4 -2 -2

>[4,] 0 0 0 0

> > b

> [,1]

>[1,] 0

>[2,] 2

>[3,] 2

>[4,] 0

>

>(This is ex Ch1, 2.2 of Artin, Algebra).

>So, 3 eqs in 4 unknowns. One can easily use row-reductions to find a

>homogeneous solution(b=0) of:

>X_1 = 0, X_2 = -c/2, X_3 = c, X_4 = 0

>

>and solutions of the above system are:

>X_1 = 2/3, X_2 = -1/3-c/2, X_3 = c, X_4 = 0.

>

>So the Kernel is 1-D spanned by X_2 = -X_3 /2, (nulliity=1), rank is 3.

>

>In R I use solve():

> > solve(a,b)

>Error in solve.default(a, b) :

> Lapack routine dgesv: system is exactly singular

>

>and it gives the error that the system is exactly singular, since it seems

>to be trying to invert a.

>So my question is:

>Can R only solve non-singular linear systems? If not, what routine should I

>be using? If so, why? It seems that it would be simple and useful enough to

>have a routine which, given a system as above, returns the null-space

>(kernel) and the particular solution.

>

>

>

>

>--

>View this message in context:

>

http://r.789695.n4.nabble.com/Finding-solution-set-of-system-of-linear-equations-tp3541490p3541490.html>Sent from the R help mailing list archive at Nabble.com.

>

>______________________________________________

>

[hidden email] mailing list

>

https://stat.ethz.ch/mailman/listinfo/r-help>PLEASE do read the posting guide

http://www.R-project.org/posting-guide.html>and provide commented, minimal, self-contained, reproducible code.

================================================================

Robert A. LaBudde, PhD, PAS, Dpl. ACAFS e-mail:

[hidden email]
Least Cost Formulations, Ltd. URL:

http://lcfltd.com/824 Timberlake Drive Tel: 757-467-0954

Virginia Beach, VA 23464-3239 Fax: 757-467-2947

"Vere scire est per causas scire"

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide

http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.