Ex17-12-20, R - rextester

Language: Layout:

# Gráfico em 3 dimensões

# Curva
# Superfície z=f(x,y)

x <- seq(0,2.5, by=0.05)
y <- seq(0,1,by=0.05)
mf<-function(s,r){z=sqrt(s*r)*exp(-sqrt(s^3)/2-r^2)# f(x,y); auxiliar para gráfico
z
}

require(grDevices) # for trans3d
z <- outer(x, y, mf)
z[is.na(z)] <- 4
op <- par(bg = "white")
persp(x, y, z, theta = 60, phi = 20, expand = 0.5)

# Curva (u(t),v(t))

t=seq(0,2.5,by=0.01)
u<-function(s){s*cos(s)+2}
v<-function(s){(-s^2)*sin(s)/4+1}

curve(u,0,2.5)

curve(v,0,2.5)

plot(u(t),v(t),'l',col="red")   # curva (u(t),v(t))

persp(x, y, z, theta = 60, phi = 20, expand = 0.5)->res
lines(trans3d(u(t),v(t),mf(u(t),v(t)), res), col = "red", lwd = 2) # Inclui curvas no gráfico em 3d

# Método de Euler oara resolver a equação e'(t)=grad f(e(t)), e(0)=(0.5,1)

gradmf<-function(u){ # gradiente de f(x,y)
s=u[1]
r=u[2]
dx=exp(-sqrt(s^3)/2-r^2)*(r/(2*sqrt(s*r))-3*s*sqrt(r)/4)
dy=exp(-sqrt(s^3)/2-r^2)*(s/(2*sqrt(s*r))-2*sqrt(s*r)*r)
c(dx,dy)
}

t0=0  # tempo inicial
tf=5  # t final 
e0=c(0.5,1) # condição inicial 
n=50
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áximo"); Y[,length(tt)]
print("Aproximação para o valor máximo"); 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(t,mf(u(t),v(t)),'l',col="red")   # curva f(u(t),v(t))

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

# Gradiente igual a zero

s<-seq(0,2.5,length=1000)
r<-seq(0,1,length=1000)
z<-outer(s,r,function(s,r) exp(-sqrt(s^3)/2-r^2)*(r/(2*sqrt(s*r))-3*s*sqrt(r)/4))
contour(s,r,z,levels=0) # dx=0
par(new=TRUE)
z<-outer(s,r,function(s,r) exp(-sqrt(s^3)/2-r^2)*(s/(2*sqrt(s*r))-2*sqrt(s*r)*r))
contour(s,r,z,levels=0) # dy=0

[ + ] Show input
Absolute running time: 2.44 sec, cpu time: 1.74 sec, memory peak: 105 Mb, absolute service time: 2,64 sec