ExDecaiEulerMin28-01-21, R

Language: Layout:

# 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
x=u[1]
y=u[2]
f1=y*x^3+y^4+2*x*y+x+y+1 
f1 
}

print("Teste de f1:"); f1(c(1,2))

f2<-function(u){   # Definição de f1
x=u[1]
y=u[2]
f2=y*x^2+2*x*y+x+5
f2
}

print("Teste de f2:"); f2(c(1,2))

f<-function(u){   # Definição de f=(f1,f2)
c(f1(u),f2(u))
}

print("Teste de f:"); f(c(1,2))

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

print("Teste de mf=|f|^2:"); mf(1,2)

gradmf<-function(u){ # gradiente de mf(x,y) aproximado
s=u[1]
r=u[2]
h=10^(-5)
dnfx= mf(s+h,r)-mf(s-h,r)
dnfy=mf(s,r+h)-mf(s,r-h)
p=c(dnfx,dnfy)/(2*h)
return(p)
}

print("Teste de gradmf:"); gradmf(c(1,2))

Jacf<-function(u){ # Jacobiano aproximado de f(x,y)
h=10^(-5)

v=u; v[1]=u[1]+h; v1=u; v1[1]=u[1]-h; 
dfx= f(v)-f(v1)

v=u; v[2]=u[2]+h; v1=u; v1[2]=u[2]-h; 
dfy= f(v)-f(v1)

A=matrix(0,2,2)
A[1,]=dfx; A[2,]=dfy
p=A/(2*h)
return(p)
}

print("Teste do Jacobiano de f:"); Jacf(c(1,2))

iJacf<-function(u){ # inversa do Jacobiano  aproximado de f(x,y) (caso 2x2)
A=Jacf(u); B=0*A
d=A[1,1]*A[2,2]-A[1,2]*A[2,1] # determinante
B[1,1]=A[2,2]; B[1,2]=-A[2,1]; B[2,1]=-A[1,2]; B[2,2]=A[1,1]
p=t(B)/d
return(p)
}

print("Teste da inversa do Jacobiano de f:"); iJacf(c(1,2))

iJacf(c(1,2))%*%Jacf(c(1,2))

alpha<-function(u){min(eigen(t(Jacf(u))%*%Jacf(u))$values)}

#-------------------------------------------------------------

require(grDevices) # for trans3d (colocar curva na superfície)
x <- seq(-5,2, by=0.05)    # Malha para gráficos
y <- seq(-2,2,by=0.05)
z <- outer(x, y, mf)
z[is.na(z)] <- 4
op <- par(bg = "white")
persp(x, y, z, theta = 60, phi = 20, expand = 1) # gráfico de mf
contour(x,y,z,levels=c(0.5,1,5,10))              # curvas de nível de mf

#
persp(x, y, z, theta = 60, phi = 20, expand = 1)->res # for trans3d (colocar curva na superfície)

# Método de Euler para resolver a equação u'(t)=-grad mf(u(t)), e(0)=(2,-2)

t0=0  # tempo inicial
tf=1  # t final 
e0=c(2,-2) # condição inicial 
mf(e0[1],e0[2]) # teste da escolha
n=5000
h=(tf-t0)/n  # Tamanho do passo

tt=seq(t0,tf,by=h)

Y=matrix(0,2,length(tt))
Y[,1]=e0

for ( i in 1:(length(tt)-1)){
Y[,i+1]=Y[,i]-h*gradmf(Y[,i])
}

print("Aproximação para o ponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[1,length(tt)],Y[2,length(tt)])

lines(trans3d(Y[1,],Y[2,],mf(Y[1,],Y[2,]), res), col = "blue", lwd = 2) # Inclui curvas

plot(Y[1,],Y[2,],col="blue")

plot(tt,mf(Y[1,],Y[2,]), col = "blue")

#------------------------------

points(tt,mf(Y[1,1],Y[2,1])*exp(-sqrt(alpha1)*tt), col = "red",xlim=c(0,1))

plot(tt,sqrt(alpha0), col = "green")

# --------------------------- Euler

t0=0  # tempo inicial
tf=1  # t final 
e0=c(2,-2) # condição inicial 
mf(e0[1],e0[2]) # teste da escolha
n=5000
h=(tf-t0)/n  # Tamanho do passo

tt=seq(t0,tf,by=h)

Y=matrix(0,2,length(tt))
Y[,1]=e0

for ( i in 1:(length(tt)-1)){
Y[,i+1]=Y[,i]-h*t(Jacf(Y[,i]))%*%f(Y[1,i],Y[2,i])
}

print("Aproximação para o ponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[1,length(tt)],Y[2,length(tt)])

plot(Y[1,],Y[2,],col="blue")

plot(tt,mf(Y[1,],Y[2,]), col = "blue")

alpha0=0*tt
for ( i in 1:length(tt)){   # Para teste de decaimento
alpha0[i]=alpha(Y[,i])}

points(tt,mf(Y[1,1],Y[2,1])*exp(-sqrt(alpha1)*tt), col = "red",xlim=c(0,1))

# ---------------------------Euler---->Newton

t0=0  # tempo inicial
tf=10  # t final 
e0=c(2,-2) # condição inicial 
mf(e0[1],e0[2]) # teste da escolha
n=10
h=(tf-t0)/n  # Tamanho do passo

tt=seq(t0,tf,by=h)

Y=matrix(0,2,length(tt))
Y[,1]=e0

for ( i in 1:(length(tt)-1)){
Y[,i+1]=Y[,i]-h*iJacf(Y[,i])%*%f(Y[1,i],Y[2,i])
}

print("Aproximação para o ponto de mínimo"); Y[,length(tt)]
print("Aproximação para o valor mínimo"); mf(Y[1,length(tt)],Y[2,length(tt)])

plot(Y[1,],Y[2,],col="blue")

plot(tt,mf(Y[1,],Y[2,]), col = "blue")

[ + ] Show input
Absolute running time: 0.49 sec, cpu time: 0.76 sec, memory peak: 48 Mb, absolute service time: 0,7 sec