3

Attention:

The question was modified for numerical solution. (see Update 2)

I have the following system of differential equations:

\begin{align} \dot x & = -x+F(a\,y+b),\\ \dot y & = -y +F(a\,x+c),\\ \text{for} \quad F(x) & = \dfrac{1}{1+e^{-x}}. \end{align}

I need to draw the landscape like it is done for gravitation:

enter image description here

I am rather far from physics, so I have two questions:

1) Is it possible to calculate potential analytically?

In case it is impossible:

2) Could You explain me (in a simple way) how I may calculate potential numerically if I know only dx and dy?

Note: If it will help, $x$ and $y$ are in the range from 0 to 1, and potential in point (0,0) is some constant, for example, zero.

I understand, that it sounds like a homework, but it is not, I would be really thankful for an explanation.

Update1:

Ok, I understand, that there is no analytical solution. But If I will do the following:

I may calculate numerically for a given point $(x_0,y_0)$ and very small $\Delta t$,

$ x_1 = x(x_0,y_0,t_0+\Delta t),x_2 = x(x_0,y_0,t_0+2\Delta t),x_3 = x(x_0,y_0,t_0+3\Delta t)$

and the same for y:

$ y(x_0,y_0,t_0+\Delta t),y(x_0,y_0,t_0+2\Delta t),y(x_0,y_0,t_0+3\Delta t)$,

then if I understand correctly I may calculate acceleration as:

$$a_x = ((x_3-x_2)/\Delta t)-(x_2-x_1)/\Delta t))/\Delta t$$

$$a_y = ((y_3-y_2)/\Delta t)-(y_2-y_1)/\Delta t))/\Delta t$$

Is it possible to draw a landscape in this case? I do not take care of precise values, so I am indifferent in scaling. I just need hills and valleys for this.

Update 2:

So, thank you very much for Your answers, but they did not work in my hands at least. Here is my attempt to solve my problem numerically (maybe it will became clear what I want, and why other approaches did not work for me).

So I may code my system in the following way (Matlab):

a = -20;
b = 10;
c = 10;

% define functions
F = @(x) 1./(1+exp(-x));
delta_x = @(x,y) -x+F(a.y+b);
delta_y = @(x,y) -y+F(a.x+c);


% get grid
nBins = 20;
x_range = linspace(-0.5,1.5,nBins);
y_range = linspace(-0.5,1.5,nBins);
[x,y]=meshgrid(x_range,y_range);


% get dx,dy
dx =  delta_x(x,y);
dy =  delta_y(x,y);

Let's plot vector field (I publish picture for smaller bin number):

figure()
quiver(x,y,dx,dy)
xlabel('x')
ylabel('y')

enter image description here

You may see, that I have a symmetric field (no surprise with symmetric parameters), where I have two stable wells at (1,0) and (0,1), two hills upper-right and lower left corners, and an unstable saddle in (0.5,0.5). Now I just want to visualise this landscape in 3D. The first (even I understand, that incorrect) approach, just to plot absolute values of velocities (speed):

Z = sqrt(dx.^2+dy.^2);
surf(x,y,Z);

enter image description here

Almost perfect picture, but at the point (0.5,0.5) I have stable well, instead of unstable saddle.

Now I try to integrate the field:

Z = cumsum(dx,2)+cumsum(dy,1);

The result is ugly: it has no symmetry, no wells...

enter image description here

Could You help me to find a mistake?

Hope, I just integrate in a wrong way...

Community
  • 1
zlon
  • 131

4 Answers4

3

Drawing a potential landscape "like it is done for gravitation" requires a second-order differential equation, of the form \begin{align} \frac{\mathrm d^2}{\mathrm dt^2}x(t) & = f(x(t),y(t)) \\ \frac{\mathrm d^2}{\mathrm dt^2}y(t) & = g(x(t),y(t)). \end{align} Your system of equations does not have this form, so the concept of a potential landscape does not make sense for your system.

Emilio Pisanty
  • 132,859
  • 33
  • 351
  • 666
  • And it is not possible to write equations in this form? – zlon Dec 12 '18 at 16:56
  • @zion Apologies, I posted without really realizing what I had written. See edited answer. – Emilio Pisanty Dec 12 '18 at 17:17
  • In other words to calculate the potential I need to know the force, but I have only velocity. If it is correct, why it is not possible to differentiate velocity to obtain a force? – zlon Dec 12 '18 at 18:09
  • 4
    Because if you have a system of the form $\dot x = g(x)$ and you differentiate it, you will get an acceleration of the form $\ddot x = \dot x g'(x)$, which depends on the velocity, and that cannot be shaped in the form $g(x) = \nabla U(x)$, i.e. forces that come from potentials cannot depend on the velocity. – Emilio Pisanty Dec 12 '18 at 21:17
1

I'm not sure it is appropriate to be solving your ODEs for you, but the system \begin{align} d x/dt & = -x+F(a\,y+b),\\ dy/dt & = -y +F(a\,x+c) \end{align} may have t decoupled by dividing the two equations. That is, for all t, you must have $$ \frac{dy}{dx}=\frac{-y+F(ax+c)}{-x+F(ay+b)}~~, $$ whose integration might, or might not, be straightforward, to yield $y(x)$, the same function for each t.

I cannot see an obvious symmetry (yet) between x and y, given your differing constants b,c.

Your parameter a is superfluous, and could be absorbed into x,y, so $$ \frac{dy}{dx}=\frac{-y+aF(x+c)}{-x+aF(y+b)}~~. $$

On the other hand, it might be conceivable you have no interest in the solutions, and you are only interested in plotting the potential leading up to them, in which case you simply plot, for constant E, $$ V(x,y)=E -\frac{m}{2}\left ((F(ay+b)-x)^2+(F(ax+c)-y)^2 \right ). $$ Again, lacking axial symmetry. Would this be what you might be after?

Cosmas Zachos
  • 62,595
1

It is also possible to look for a "landscape" in which the system always moves "downhill". In mathematical terms, supposed we have $\mathbf{v} = (x,y)$, and the equations of motion are $\dot{\mathbf{v}} = \mathbf{f}$. If we can find a function $U(x,y)$ such that $\mathbf{f} = -\nabla U$, then the system will always move in the direction that $U$ is decreasing, reaching equilibrium when $U$ is at a local extremum. This may correspond to your intuitive notion of "how it is done for gravity" (though as Emilio points out, the situation's a bit more nuanced than that.)

However, for such a $U$ to exist, we have to have $\nabla \times \mathbf{f} = 0$ for all values of $x$ and $y$, since the curl of a gradient is always zero. In your case, you have $$ \mathbf{f} = \begin{bmatrix} - x + F(ay + b) \\ -y + F(ax + c) \end{bmatrix} $$ and the curl of this quantity is $$ (\nabla \times \mathbf{f})_z = \frac{\partial f_y}{\partial x} - \frac{\partial f_x}{\partial y} = a \left( F'(ax + c) - F'(ay + b) \right) $$ Since this must be true for all values of $x$ and $y$, it must be the case that $F'(z)$ is a constant for all values of its argument $z$, which implies that $F(z) = \alpha z + \beta$ for some values of $\alpha$ and $\beta$. For any other form of the function $F$ (including the form given in the OP), there will be some point $(x,y)$ at which $\nabla \times \mathbf{f} \neq 0$, and thus we cannot find a potential $U$ with the desired property.

EDIT: Note that in the revised question, the OP is effectively calculating and plotting $|\mathbf{f}|$ in the first surface graph. In the second, they are effectively calculating $$ g(x,y) = \int_0^x f_x(x', y) dx' + \int_0^y f_y(x, y') dy'. $$ If the potential $U$ exists, it could be found via a similar technique, but the integrals would have to be subtly different; you'd have to have $$ U(x,y) \overset{?}{=} \int_0^x f_x(x', 0) dx' + \int_0^y f_y(x, y') dy' $$ or $$ U(x,y) \overset{?}{=} \int_0^x f_x(x', y) dx' + \int_0^y f_y(0, y') dy'. $$ These two quantities correspond to taking $\int_{(0,0)}^{(x,y)} \mathbf{f} \cdot d\mathbf{r}$ via the points $(0,x)$ and $(0,y)$, respectively. But for $U$ to be well-defined, these two quantities to be equal; and for that to happen, you have to have $\nabla \times \mathbf{f} = 0$, which you don't have in this case.

Michael Seifert
  • 48,772
1

FWIW, it is in principle possible to (locally) write OP's system of two coupled first order ODEs as a Hamiltonian system

$$ \dot{x}~=~\{x,H\}, \qquad \dot{y}~=~\{y,H\}, $$

with a (non-canonical) symplectic structure $\{\cdot,\cdot\}$ and a Hamiltonian $H(x,y)$, cf. this Phys.SE post. The integral curves are then constant-energy contours $H(x,y)=E$ in a potential landscape $(x,y)\mapsto H(x,y)$.

Qmechanic
  • 201,751