Exercicio3

Joao_Bianco · updated July 15, 2021 · R
# Gráfico em 3 dimensões

# Resolver a equação f(x,y)=(f1(x,y),f2(x,y))=(0,0).

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

x <- seq(-4,4, by=0.05)
y <- seq(-4,4,by=0.05)
mf1<-function(r,s){u=c(r,s); f1(u)}# auxiliar para gráfico
z=matrix(0,length(x),length(y))
for ( i in 1:length(x)){
  for ( j in 1:length(y)){z[i,j]=f1(c(x[i],y[j]))}}


require(grDevices) # for trans3d
persp(x, y, z, theta = 30, phi = 40, expand = 1)
contour(x,y,z,level=0,col="red")
par(new=TRUE)

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

zz=matrix(0,length(x),length(y))
for ( i in 1:length(x)){
  for ( j in 1:length(y)){zz[i,j]=f2(c(x[i],y[j]))}}

contour(x,y,zz,level=0,col="blue")

legend(2,3.5, legend=c("f1(u)", "f2(u)"),
       col=c("red","blue"), lty=c(1,1), cex=0.8)

persp(x, y, zz, theta = 30, phi = 40, expand = 1)

f<-function(u){c(f1(u),f2(u))} # f(x,y)

mf<-function(u){   # Definição de mf=|f|^2; com coordenadas para gráfico
  f1(u)^2+f2(u)^2 
}


#---------#

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

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

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

Jf<-function(u){
  x=u[1];y=u[2]
  p=matrix(0,2,2)
  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)
  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 )
  p[1,]=df1
  p[2,]=df2
  p
}

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

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

ZeroEuler<-function(u0,t,n){
  u=matrix(0,2,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 

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

x=u[1,]; y=u[2,]
plot(x,y,'l',col="blue", main="Euler-Newton para Jf(u(t))u'(t)=-alpha f(u(t)), com n=10")

Comments

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