how to calculate this simultaneous equations in r

Question

I am trying to solve the following two equations simultaneously for a and b in R:

0.1 = a /(1+a) * 0.9639 + a(1+b) / (1+a(1+b)) * 0.0324 + a(1+b)^2 / (1+a(1+b)^2) * 0.0036 + a(1+b)^4 / (1+a(1+b)^4) * 0.0001

0.03 = [(a/(1+a)-0.1)^2 * 0.9639 + (a(1+b) / (1+a(1+b))-0.1)^2 * 0.0324 +(a(1+b)^2 / (1+a(1+b)^2)-0.1)^2 * 0.0036 + (a(1+b)^4/(1+a(1+b)^4)-0.1)^2*0.0001]/0.09

I have tried to use the rootSolve package to find a solution and got an error message.

The code that I used was as follows:

library(rootSolve)

model=function(x){ 
    f1=((x[1]/(1+x[1]))*0.9639+((x[1]*(1+x[2]))/(1+x[1]*(1+x[2])))*0.0324+((x[1]*(1+x[2])^2)/(1+x[1]*(1+x[2])^2))*0.0036+(x[1]*(1+x[2])^4)/(1+x[1]*(1+x[2])^4)*0.0001)-0.1
    f2=(((x[1]/(1+x[1])-0.1)^2*0.9639+(x[1]*(1+x[2])/(1+x[1]*(1+x[2]))-0.1)^2*0.0324+(x[1]*(1+x[2])^2/(1+x[1]*(1+x[2])^2)-0.1)^2*0.0036+(x[1]*(1+x[2])^4/(1+x[1]*(1+x[2])^4)-0.1)^2*0.0001)/0.09)-0.03
    c(f1=f1,f2=f2)
} 

solution=multiroot(f=model, start=c(1,1))

Please be more clear, with your problem. Also mention hat you have tried ? — vrajs5, May 27 '14 at 07:04
Welcome to StackOverflow! The fist thing to get high quality answers is to write high quality questions around here. You might want to read http://stackoverflow.com/questions/how-to-ask. Apart from that you might find some hints for your problem here: http://stackoverflow.com/questions/18958372/solve-systems-of-non-linear-equations-in-r-black-scholes-merton-model — Thilo, May 27 '14 at 07:04
R isn't really suited for algebra (although there may be a package for it, or possibly you could use an optimiser function). Try maxima: http://maxima.sourceforge.net/ — Spacedman, May 27 '14 at 07:08
Now include the error message you get when you try that. Your problem is (first) that you haven't learnt how to write expressions in R properly. `x[1](1+x[2])` is **NOT** a multiplication!!! — Spacedman, May 27 '14 at 07:44
I'm sorry for my mistakes~ . I've modified it , and it gave an answer but it is a wrong one. I wonder if my starting value is not correct. — user3678637, May 27 '14 at 08:28

score 2 · Accepted Answer · answered May 27 '14 at 07:57

After fixing (and slightly tidying up) your model function:

model=function(x){
    f1=((x[1]/(1+x[1]))*0.9639+
        ((x[1]*(1+x[2]))/(1+x[1]*(1+x[2])))*0.0324 +
        ((x[1]*(1+x[2])^2)/(1+x[1]*(1+x[2])^2))*0.0036+
        (x[1]*(1+x[2])^4)/(1+x[1]*(1+x[2])^4)*0.0001)-0.1
    f2=(((x[1]/(1+x[1])-0.1)^2*0.9639 +
         (x[1]*(1+x[2])/(1+x[1]*(1+x[2]))-0.1)^2*0.0324 +
         (x[1]*(1+x[2])^2/(1+x[1]*(1+x[2])^2)-0.1)^2*0.0036+
         (x[1]*(1+x[2])^4/(1+x[1]*(1+x[2])^4)-0.1)^2*0.0001)/0.09)-0.03
    c(f1=f1,f2=f2)
}

I can run multiroot without error:

> multiroot(model,start=c(1,1))
diagonal element is zero 
[1] 2
$root
[1] -1.827359e+00  1.749761e+06

$f.root
       f1        f2 
 2.065032 47.916572 

$iter
[1] 3

$estim.precis
[1] 24.9908

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

BUT the warnings tell me it probably didn't find a true solution. For these problems you do have to experiment with starting values, especially since your problem is VERY non linear. After a bit of experimenting with calling model(c(a,b)) for various values of a and b to see which direction I might find zeroes, I hit on this:

> multiroot(model,start=c(.1,3.5))
$root
[1] 0.09989584 3.44808411

$f.root
           f1            f2 
-6.647460e-13  3.676781e-13 

$iter
[1] 3

$estim.precis
[1] 5.162121e-13

Which gives no warnings and claims to be pretty close to zero. Your answer is in the $root part.

There's possibly other solutions - you have x[2] up to powers of 4, so that could mean four solutions (some of which may be complex numbers).

If you can get an algebra system to solve the equations exactly (replace your decimal numbers with constants) then you might get a closed-form for the solutions.

Can you click the tick mark to show my answer answered your question? — Spacedman, May 27 '14 at 09:29

how to calculate this simultaneous equations in r

1 Answers1