Atividade extra 2 - Exercicio 2

Joao_Bianco · updated August 04, 2021 · R
# --------------------------------------------------------------------------------------------------
#  ------------------ 2. Resolva o Exercício 2.1.22 das notas de Aula.------------------------------
# --------------------------------------------------------------------------------------------------

f1<-function(u){   # Definição de f1
    x=u[1]
    y=u[2]
    z=u[3]
    return(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)
}

f2<-function(u){   # Definição de f2
    x=u[1]
    y=u[2]
    z=u[3]
    return(6*y^5-8*y^3+2*y-12*x*y+12*x^2*y^3+6*x^4*y+z)
}

f3<-function(u){    # Definição de f3
    x=u[1]
    y=u[2]
    z=u[3]
    return(4*z^3-8*z+x+y)
}

#--- Minimo de |grad f(u)|^2 para resolver equação f(u)=0.

f<-function(u){ # Definição de função f(u), para resolver a equação f(u)=0, com f(u)=grad gf(u).
    c( f1(u), f2(u),f3(u))
}

g<-function(u){
    t(f(u)) %*% f(u)/2
} 

Jf<-function(u){ # A matriz do jacobiano de f(u) é a matriz hessiana de gf(u).
    x=u[1]
    y=u[2]
    z=u[3]
    p=matrix(0,3,3)
    df1=c(6*y^4+2+12*x+36*x^2*y^2+24*x^2+30*x^4, -12*y+24*x*y^3+24*x^3*y,0)
    df2=c(-12*y+24*x*y^3+24*x^3*y, 30*y^4-24*y^2+2-12*x+36*x^2*y^2+6*x^4,0 )
    df3=c(1,1,12*z^2-8)
    p[1,]=df1
    p[2,]=df2
    p[3,]=df3
    p
}

gradg<-function(u){  
    p=t(Jf(u))%*%f(u)
    t(p)
}


# Método de Euler-Newton para a Equação Jf(u(t))u'(t)=-alpha f(u(t))
alpha=9.8

ZeroEuler<-function(u0,t,n){
    u=matrix(0,3,n+1)
    u[,1]=u0
    tj=t/n
    for (i in 1:n){
        h=solve( Jf(u[,i]),-tj*alpha*f(u[,i]))    # Resolve o sistema linear Jf(u)h=-tj alpha f(u)
        u[,i+1]=u[,i]+h
    }
    u
}

# Teste 1

print("Resolucao equacao 1")
u0=c(0,0,0); t=1; n = 10
u=ZeroEuler(u0,t,n) ; u [,n+1]     # Aproximadamente u(1)
f(u[,n+1])
g(u[,n+1])

A=Jf(u[,n+1]);A
A[1:2,1:2]

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

# O cálculo dos coeficientes do polinômio é feito através  das linhas de comando
	
coef<-function(A){ 
    n=length(A[1,]);	a=0*1:(n+1)
    s=0*1:n; s[1]=tr(A)
    a[1]=1;a[2]=-s[1]; B=A
    for ( i in 2:n){
        B=B%*%A; s[i]=tr(B); q=0
        for (j in 1:i){
            q=q+a[j]*s[i+1-j]
        }; a[i+1]=-q/i
    }
    a=(-1)^n*a;a
}
a=coef(A);a	

# Nesse caso  temos policarac(x)=a[1]x^3 + a[2]x^2+ a[3]x+a[4].
policarac<-function(s){
    a[1]*s^3 + a[2]*s^2+a[3]*s+a[4]
}
curve(policarac,-8,17)

# Teste 2

print("--------------------------------------------------------------------------------------------------")
print("Resolucao equacao 2")
u0=c(1,2,3); t=1; n = 10
u=ZeroEuler(u0,t,n) ; u [,n+1]     # Aproximadamente u(1)
f(u[,n+1])
g(u[,n+1])

A=Jf(u[,n+1]);A
A[1:2,1:2]

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

# O cálculo dos coeficientes do polinômio é feito através  das linhas de comando
	
coef<-function(A){ 
    n=length(A[1,]);	a=0*1:(n+1)
    s=0*1:n; s[1]=tr(A)
    a[1]=1;a[2]=-s[1]; B=A
    for ( i in 2:n){
        B=B%*%A; s[i]=tr(B); q=0
        for (j in 1:i){
            q=q+a[j]*s[i+1-j]
        }; a[i+1]=-q/i
    }
    a=(-1)^n*a;a
}
a=coef(A);a	

# Nesse caso  temos policarac(x)=a[1]x^3 + a[2]x^2+ a[3]x+a[4].
policarac<-function(s){
    a[1]*s^3 + a[2]*s^2+a[3]*s+a[4]
}
curve(policarac,-8,20)

Comments

Please sign up or log in to contribute to the discussion.