programación de flujo de gradiente (matlab)

Question

Considere $$f(x,y)= x \ y \ e^{-x^2-y^2}$$

Necesito encontrar el flujo gardiente de f empezando en el punto (0,1) y luego (2,4). Sé que esto implica resolver la siguiente EDO $$\dot X(t)= - \nabla f(X(t))$$

Estoy usando este código, pero estoy obteniendo resultados estresantes:

x_old = [0 0]
x_new = [2 6] 
eps = 0.01
precision = 0.000001;
fx = @(x,y) ((1-2*x.*x).*y*exp(-x.^2 - y.^2));
fy=  @(x,y) ((1-2*y.*y).*x*exp(-x.^2 - y.^2));
while abs(x_new - x_old) > precision
    x_old = x_new;
    x_new = x_old - eps * [fx(x_old(1),x_old(2)) fy(x_old(1),x_old(2))];
end
x_new

EDITAR:

¡Ahora estoy usando ODE45 pero estoy obteniendo los datos iniciales como una solución sin importar donde empiezo!

function Fv=funsys(t,Y)
Fv(1,1)=(1-2*Y(1).*Y(1)).*Y(2).*exp(-Y(1).^2 - Y(2).^2);
Fv(2,1)=(1-2*Y(2).*Y(2)).*Y(1).*exp(-Y(1).^2 - Y(2).^2);

Entonces llama a ode45

[tv,Yv]=ode45('funsys',[0 1],[6,2], options)

¿Qué me estoy perdiendo?

Answer 1

En primer lugar, podemos realizar un análisis cualitativo del campo vectorial mediante un gráfico de campo vectorial .

Se puede observar que hay puntos fijos de repulsión cerca de $(0.6,0.6)$ y $(-0.6,-0.6)$ y los puntos fijos de atracción $(0.6,-0.6)$ y $(-0.6,0.6)$ . También hay un punto de ensillamiento en $(0,0)$ . Este gráfico se ha realizado en MATLAB con el código

[X,Y] = meshgrid(linspace(-2,2,21)); % create mesh
dfx = @(x,y) -(y-2*x.^2.*y).*exp(-(x.^2+y.^2)); % x partial derivative
dfy = @(x,y) -(x-2*x.*y.^2).*exp(-(x.^2+y.^2)); % y partial derivative
Fx = dfx(X,Y);Fy = dfy(X,Y); % calculate the partial derivatives
quiver(X,Y,Fx,Fy) % plot the vector field
axis square,axis([-1,1,-1,1]*2),xlabel('x'),ylabel('y') % set plot options

A partir de esto, esperamos que la condición inicial $(0,1)$ para tender hacia el punto fijo cerca de $(-0.6,0.6)$ . Aunque no se muestra en el gráfico, esperamos que la condición inicial $(2,4)$ para aumentar sin límites.

El siguiente código se utiliza para producir estos resultados. El sistema ODE se resuelve con ODE45. Esto está diseñado para ser ejecutado después del fragmento de código anterior.

g = @(t,y) [dfx(y(1),y(2));dfy(y(1),y(2))];
ic = [0;1];
tend = 10;
[t,sol] = ode45(g,[0,tend],ic);

% generate graph
subplot(2,1,1)
quiver(X,Y,Fx,Fy),hold on
plot(sol(:,1),sol(:,2),'r','Linewidth',2)
axis([-1,0.5,0.5,1.1])
xlabel('x'),ylabel('y')
subplot(2,1,2)
plot(t,sol(:,1),'b',t,sol(:,2),'r')
legend('x(t)','y(t)','Location','SouthEast')

El gráfico superior muestra el $(x,y)$ posición superpuesta en el campo vectorial. Los componentes individuales se muestran en función del tiempo en el gráfico inferior.

Esto se puede repetir con la segunda condición inicial $(2,4)$ . tend se ha incrementado a $10^{12}$ en este gráfico. Obsérvese también la escala logarítmica en el $x$ -eje.

La resolución del sistema ODE para ambas condiciones iniciales, como se ha hecho, confirma el análisis cualitativo realizado inicialmente.