Otavio-fez, R - rextester

Language: Layout:

#R version 3.3.2 
  
print("Hello, world!")
#R version 3.3.2 
f<-function(s){
  x=s[1];y=s[2]; z=s[3]
  
  p = y^6-2*y^4+y^2+1+2*x-6*x*y^2+3*x^2*y^4+x^2+2*x^3+3*x^4*y^2+2*x^4+x^6+z^4-4*z^2+4 + x*z + y*z
  p
}

#derivadas em x,y,z

dfx<-function(s){
  x=s[1];y=s[2]; z=s[3]
  
  p = 2 - 6*y^2 + 6*x*y^4 + 2*x + 6*x^2 + 12*x^3*y^2 + 8*x^3 + 6*x^5 + z
  p
}

dfy<-function(s){
  x=s[1];y=s[2]; z=s[3]
  
  p = 6*y^5 - 8*y^3 + 2*y - 12*x*y + 12*x^2*y^3 + 6*x^4*y + z
  p
}

dfz<-function(s){
  x=s[1];y=s[2]; z=s[3]
  
  p = 4*z^3 - 8*z + x + y
  p
}

#gradiente de f(x,y,z)

grad<-function(s){
  
  p = c(dfx(s), dfy(s), dfz(s))
  p
  
}

#resolvendo numericamente a equacao diferencial (euler)

u0=c(0,-2, 0)
#u0 = c(1,1,1)

a=0

b=1

n=1000000

h=(b-a)/n

u=u0

for ( i in 1:n){ u=u-h*grad(u)}
print('u')
u

print("f(u)")
f(u)

print("grad(u)")
grad(u)

# calculando a matriz hessiana de f(x,y,z)

Hf<-function(s){ 
x=s[1];y=s[2]; z=s[3]

H=matrix(0,3,3)

H[1,1]=6*y^4 + 2 + 12*x + 36*x^2*y^2 + 24*x^2 + 30*x^4

H[1,2]= -12*y + 24*x*y^3 + 24*x^3*y;

H[1,3] = 1

H[2,1] = H[1,2]

H[2,2]=30*y^4-24*y^2+2-12*x+36*x^2*y^2+6*x^4

H[2,3] = 1

H[3,3]=12*z^2-8

H
}

print("Matriz Hessiana de Hf(u)")
Hf(u)

# Calculando o polinomio caracteristico de Hf(x,y,z) pelas Identidades de Newton

tr<-function(A){p=0; n=length(A[1,]); for ( i in 1:n ){p=p+A[i,i]}; p}

coef<-function(A){ 
  n = length(A[1,]); # dimensao da matriz
  a = 0*1:(n+1);     # vetor de tamanho n + 1
  s = 0*1:n;         # vetor de tamanho n
  s[1] = tr(A)       # calculando o traco da matriz A
  a[1] = 1;          # elemento p0
  a[2] = -s[1];      # elemento p1 (traco da matriz A * -1)
  B = A;             
  for ( i in 2:n){   # repeticao para calcular o traco de H^i
    B=B%*%A;         # multiplicacao das matrizes
    s[i]=tr(B);      # calcular o traco da matriz H^i
    q=0       
    for (j in 1:i){  # somatorio para calcular bj
      q = q + a[j]*s[i+1-j]
    }; 
    a[i+1]=-q/i;
  }; 
  a=(-1)^n*a;        # invertendo o sinal de todos os b's  
  a
}

print("Polinomio caracteristico de Hf(u)")
coef(Hf(u))

[ + ] Show input