gradient flow programming (matlab)

Question

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

I need to find the gradient flow of $f$ starting at the point (0,1) then (2,4). I know that this entails solving the following ODE $$\dot X(t)= - \nabla f(X(t))$$

I am using this code, but am getting strange results:

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

EDIT:

Now am using ODE45 but am getting the initial data as a solution no matter where I start!

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);

Then call ode45

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

What am I missing??

Answer 1

Firstly, we can conduct qualitative analysis of the vector field by way of a vector field plot.

You can observe that there are repelling fixed points near $(0.6,0.6)$ and $(-0.6,-0.6)$ and attracting fixed points $(0.6,-0.6)$ and $(-0.6,0.6)$. There is also a saddle point at $(0,0)$. This graph was produced in MATLAB with the code:

[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

From this, we expect the initial condition $(0,1)$ to tend towards the fixed point near $(-0.6,0.6)$. While not shown on the graph, we expect the initial condition $(2,4)$ to increase without bound.

The following code is used to produce these results. The ODE system is solved using ODE45. This is designed to be execute after the above code snippit.

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')

The top graph shows the $(x,y)$ position overlaid on the vector field. The individual components are shown w.r.t. time in the bottom graph.

This can be repeated with the second initial condition $(2,4)$. tend has been increased to $10^{12}$ in this graph. Also note the log scale on the $x$-axis.

Solving the ODE system for both initial conditions, as has been done, confirms the qualitative analysis initially performed.