Newton-Raphson

newton-raphson

c1= 0.35; c2=0.14
old=c(0,0.20);new=U=I=NULL;
f <- function(old){
w0=old[1];w1=old[2];
f1 = c1*(w0^2+w1^2+(1-w0-w1)^2) -w1*(1-w1)
f2 = c2*(w0^2+w1^2+(1-w0-w1)^2) -w0*(1-w0-w1)
U=c(f1,f2)
d11 = c1*(2*w0-2*(1-w0-w1));
d12 = c1*(2*w1-2*(1-w0-w1))-1+2*w1;
d21 = c2*(2*w0-2*(1-w0-w1))-1+2*w0+w1;
d22 = c2*(2*w1-2*(1-w0-w1))+w0;
I=matrix(c(d11,d21,d12,d22),nrow=2)
new = old-solve(I)%*%U
return(new)
}

isBad=TRUE;count=0;
while(isBad){
new = f(old)
  if( sum((old-new)^2)<10^(-8) ){
        isBad= FALSE;
    }else{
        old=new
        count=count+1;
    }
}
new;count;

##        [,1]
## [1,] 0.1065
## [2,] 0.2294

## [1] 3

w = c(new[1],new[2],1-sum(new))

n1 = w[1]*w[2]+w[2]*w[3];
n2 = w[1]*w[3];
d = sum(w^2);

n1/d;n2/d;

## [1] 0.35

## [1] 0.14