Search code examples
rsystemnonlinear-functions

Solve system of non-linear equations


I am trying to solve the following system of four equations. I have tried using the "rootSolve" package but it does not seem like I can find a solution this way.

The code I am using is the following:

model <- function(x) {
F1 <- sqrt(x[1]^2 + x[3]^2) -1
F2 <- sqrt(x[2]^2 + x[4]^2) -1
F3 <- x[1]*x[2] + x[3]*x[4]
F4 <- -0.58*x[2] - 0.19*x[3]
c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
}
(ss <- multiroot(f = model, start = c(0,0,0,0)))

But it gives me the following error:

Warning messages:
1: In stode(y, times, func, parms = parms, ...) :
error during factorisation of matrix (dgefa);         singular matrix
2: In stode(y, times, func, parms = parms, ...) : steady-state not reached

I have changed the starting values, as suggested in another similar answer, and for some I can find a solution. However, this system - according to the source I am using - should have an uniquely identified solution. Any idea about how to solve this system?

Thanks you!


Solution

  • Your system of equations has multiple solutions. I use a different package to solve your system: nleqslv as follows:

    library(nleqslv)
    
    model <- function(x) {
       F1 <- sqrt(x[1]^2 + x[3]^2) - 1
       F2 <- sqrt(x[2]^2 + x[4]^2) - 1
       F3 <- x[1]*x[2] + x[3]*x[4]
       F4 <- -0.58*x[2] - 0.19*x[3]
       c(F1 = F1, F2 = F2, F3 = F3, F4 = F4)
    }
    
    #find solution
    xstart  <-  c(1.5, 0, 0.5, 0)
    nleqslv(xstart,model)
    

    This gets the same solution as the answer of Prem.

    Your system however has multiple solutions. Package nleqslv provides a function to search for solutions given a matrix of different starting values. You can use this

    set.seed(13)
    xstart <- matrix(runif(400,0,2),ncol=4)
    searchZeros(xstart,model)
    

    (Note: different seeds may not find all four solutions)

    You will see that there are four different solutions:

    $x
         [,1]          [,2]          [,3] [,4]
    [1,]   -1 -1.869055e-10  5.705536e-10   -1
    [2,]   -1  4.992198e-13 -1.523934e-12    1
    [3,]    1 -1.691309e-10  5.162942e-10   -1
    [4,]    1  1.791944e-09 -5.470144e-09    1
    .......
    

    This clearly suggests that the exact solutions are as given in the following matrix

    xsol <- matrix(c(1,0,0,1,
                     1,0,0,-1,
                    -1,0,0,1,
                    -1,0,0,-1),byrow=TRUE,ncol=4)
    

    And then do

    model(xsol[1,])
    model(xsol[2,])
    model(xsol[3,])
    model(xsol[4,])
    

    Confirmed! I have not tried to find these solutions analytically but you can see that if x[2] and x[3] are zero then F3 and F4 are zero. The solutions for x[1] and x[4] can then be immediately found.