数値計算の技法 星野聴先生 コロナ社

#数値計算の技法　p.61

#fiacco,mccormick(1968),murtagh,sargent(1969)

f=function(x){

f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2

return(c(f1,f2))

}

X=rep(1,2)

H=diag(1,length(X))

ite=10000

eta=10^(-2)

h=0.01

for(l in 1:ite){

f_val=f(X)

sigma=t(t(-H%*%f_val))

X=X+eta*sigma

y=f(X)-f_val

H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

print(f(X))
 
}

#数値計算の技法　p.49

#stewart(1967)

f=function(x){

f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2

return(c(f1,f2))

}

X=rep(1,2)

H=diag(1,length(X))

ite=10000

eta=10^(-2)

h=0.01

for(l in 1:ite){

f_val=f(X)

sigma=t(t(-H%*%f_val))

X=X+eta*sigma

y=f(X)-f_val

H=H+eta*sigma%*%t(sigma)/as.numeric(t(sigma)%*%y)-H%*%y%*%t(y)%*%H/as.numeric(t(y)%*%H%*%y)

print(f(X))
 
}

#数値計算の技法　p.107

#levenberg(1944)

f=function(x){

f1=x[1]^2+x[2]^2-1

f2=x[1]^2-x[2]^2

return(c(f1,f2))

}

X=rep(1,2)

ite=10000

eta=10^(-2)

h=0.01

w=0.1

for(l in 1:ite){

df=array(0,dim=c(length(X),length(X)))

for(j in 1:length(X)){
  
vec=X;vec[j]=vec[j]+h  
  
df[j,]=(f(vec)-f(X))/h  
  
}
  
A=df%*%t(df)+diag(w,length(X))  

y=df%*%f(X)

dx=svd(A)$u%*%diag(1/svd(A)$d)%*%t(svd(A)$v)%*%(-y)

X=X+dx

print(f(X))
 
}