Re: [R] How to solve a non-linear system of equations using R

From: Ravi Varadhan <RVaradhan_at_jhmi.edu>
Date: Wed, 04 Jun 2008 14:55:53 -0400

Jorge,

You can use the package "BB" to try and solve this problem.

I have re-written your functions a little bit.

 # ------------------------------

# Constants
# ------------------------------
l=1
m=0.4795
s=0.4795


# ------------------------------
# Functions to estimate f_i-k_i
# ------------------------------

myfn <- function(d){

d1 <- d[1]
d2 <- d[2]
d3 <- d[3]
d4 <- d[4]

res <- rep(NA, 4)
res[1] <-
2*d1+2*sqrt(2)*d1*d2+2*sqrt(3)*d2*d3+4*d3*d4-l*m*(1+d1^2+d2^2+d3^2+d4^2) res[2] <-
2*sqrt(2)*d2+2*d1^2+2*sqrt(6)*d1*d3+4*d2^2+4*sqrt(3)*d2*d4+6*d3^2+8*d4^2-l*( m^2+m^3*s^(-1))*(1+d1^2+d2^2+d3^2+d4^2)
res[3] <-
6*d1+12*sqrt(2)*d1*d2+18*sqrt(3)*d2*d3+48*d3*d4+2*sqrt(6)*d3+4*sqrt(6)*d1*d4 -l*(m^3+3*m^4*s^(-1)+3*m^6*s^(-2))*(1+d1^2+d2^2+d3^2+d4^2) res[4] <-
12*sqrt(2)*d2+12*d1^2+36*d2^2+72*d3^2+120*d4^2+20*sqrt(6)*d1*d3+56*sqrt(3)*d 2*d4+4*sqrt(6)*d4-l*((m^4+6*m^6*s^(-1)+15*m^6*s^(-2)+15*m^7*s^(-3))-3*l^(2)* (m^2+m^3*s^(-1))^2)*(1+d1^2+d2^2+d3^2+d4^2) res
}

myfn.opt <- function(d){
# re-writing "myfn" to be used for minimization

d1 <- d[1]
d2 <- d[2]
d3 <- d[3]
d4 <- d[4]

res <- rep(NA, 4)
res[1] <-
2*d1+2*sqrt(2)*d1*d2+2*sqrt(3)*d2*d3+4*d3*d4-l*m*(1+d1^2+d2^2+d3^2+d4^2) res[2] <-
2*sqrt(2)*d2+2*d1^2+2*sqrt(6)*d1*d3+4*d2^2+4*sqrt(3)*d2*d4+6*d3^2+8*d4^2-l*( m^2+m^3*s^(-1))*(1+d1^2+d2^2+d3^2+d4^2)
res[3] <-
6*d1+12*sqrt(2)*d1*d2+18*sqrt(3)*d2*d3+48*d3*d4+2*sqrt(6)*d3+4*sqrt(6)*d1*d4 -l*(m^3+3*m^4*s^(-1)+3*m^6*s^(-2))*(1+d1^2+d2^2+d3^2+d4^2) res[4] <-
12*sqrt(2)*d2+12*d1^2+36*d2^2+72*d3^2+120*d4^2+20*sqrt(6)*d1*d3+56*sqrt(3)*d 2*d4+4*sqrt(6)*d4-l*((m^4+6*m^6*s^(-1)+15*m^6*s^(-2)+15*m^7*s^(-3))-3*l^(2)* (m^2+m^3*s^(-1))^2)*(1+d1^2+d2^2+d3^2+d4^2) sum(res^2)
}

library(BB)  

p0 <- runif(4, -1,0)
ans1 <- dfsane(par=p0, fn=myfn)
ans2 <- spg(par=p0, fn=myfn.opt)
ans1
ans2

Note that the above does not produce a redual of zero, so the system can't be solved exactly. I tried a large number of random starting values without improving upon the solution provided by "spg". So, you may want to check your system for its correctness.

Hope this helps,
Ravi.

-----Original Message-----
From: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] On Behalf Of Jorge Ivan Velez
Sent: Tuesday, June 03, 2008 5:09 PM
To: R mailing list
Subject: [R] How to solve a non-linear system of equations using R

Dear R-list members,

I've had a hard time trying to solve a non-linear system (nls) of equations which structure for the equation i, i=1,...,4, is as follows:

                f_i(d_1,d_2,d_3,d_4)-k_i(l,m,s) = 0    (1)


In the expression above, both f_i and k_i are known functions and l, m and s are known constants. I would like to estimate the vector d=(d_1,d_2,d_3,d_4) which is solution of (1). Functions in R to estimate f_i-k_i are at the end of this message.

Any help/suggestions/comments would be greatly appreciated.

Thanks in advance,

Jorge

# ------------------------------
# Constants
# ------------------------------

l=1
m=0.4795
s=0.4795


# ------------------------------
# Functions to estimate f_i-k_i
# ------------------------------
f1=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]

res=2*d1+2*sqrt(2)*d1*d2+2*sqrt(3)*d2*d3+4*d3*d4-l*m*(1+d1^2+d2^2+d3^2+d4^2) res
}
f2=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]

res=2*sqrt(2)*d2+2*d1^2+2*sqrt(6)*d1*d3+4*d2^2+4*sqrt(3)*d2*d4+6*d3^2+8*d4^2 -l*(m^2+m^3*s^(-1))*(1+d1^2+d2^2+d3^2+d4^2) res
}
f3=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]

res=6*d1+12*sqrt(2)*d1*d2+18*sqrt(3)*d2*d3+48*d3*d4+2*sqrt(6)*d3+4*sqrt(6)*d 1*d4-l*(m^3+3*m^4*s^(-1)+3*m^6*s^(-2))*(1+d1^2+d2^2+d3^2+d4^2) res
}
f4=function(d){
d1=d[1]
d2=d[2]
d3=d[3]
d4=d[4]

res=12*sqrt(2)*d2+12*d1^2+36*d2^2+72*d3^2+120*d4^2+20*sqrt(6)*d1*d3+56*sqrt( 3)*d2*d4+4*sqrt(6)*d4-l*((m^4+6*m^6*s^(-1)+15*m^6*s^(-2)+15*m^7*s^(-3))-3*l^ (2)*(m^2+m^3*s^(-1))^2)*(1+d1^2+d2^2+d3^2+d4^2) res
}

        [[alternative HTML version deleted]]



R-help_at_r-project.org 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.

R-help_at_r-project.org 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. Received on Wed 04 Jun 2008 - 20:09:25 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 04 Jun 2008 - 20:31:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive