Good afternoon,
I have a system of $3$ highly non linear differential equations, which I have to integrate form a starting singular point $x^1=[1,1,1]$, and theoretically I have to arrive to an attractor $x^2$.
The problem comes when I try to integrate from a perturbed point $x^1 + \epsilon$, where $\epsilon$ is a random vector of norm little enough $10^{-10}$. What I get is a set of solutions which go to infinity and do not go to the point $x^2$. And if a solution ever seems clearly to go to $x^2$, it diverges at some point of the path.
Any ideas?
P.D. The system depends on a parameter $M$. When $M=1.1$, it is less stiff, and I do manage to find a point near $x^1$ which goes and arrives to the exact theoretically predicted sink point. But when $M$ is $1.2$ or more, it gets SO stiff... any ideas are welcome :) P.P.D. I have all analytical information possible: Jacobian at each point, etc. P.P.P.D: Thanks for the editing and adding the "$" :D The equations are:
$\dot{x}(1)=x(1)-1+\epsilon(x(3)-1)+\frac{1}{(\gamma M^2)}(\frac{x(2)}{x(1)}-1)$ $\dot{x}(2)=(\frac{2}{3}P(\gamma-1)M^2)(-(x(1)-1)^2+\epsilon(x(3)-2x(1)+1)(x(3)-1)+\frac{2}{\gamma(\gamma-1)M^2}(x(2)+(\gamma-1)x(1)-\gamma))$
$\dot{x}(3)=\frac{x(2)(x(1)-x(3))}{x(1)x(3)F}$
Where $P=0.67,F=17.79M^2,\gamma=\frac{5}{3},\epsilon=0.5,M\in(0,2)$ are parametres, all fixed.
As additional information, there are 2 positive eigenvalues and 1 negative at $x^1=[1,1,1]$. The ratio between the "spurious" eigenvalue and the good one (which leads to $x^2$) is: When M=1, 10. When M=1.1, 11. When M=1.2, 16. When M=1.4, 84. Etc.
The following images show how sensitive is the system to M, specifically the system integration from random points $10^{-10}$ away from $x_1$. In the first case, $M=1.2$, the right path is not found. In the second case, $M=1.1$, some points find the right path, and eventually one reaches $x^2$.
Are you sure you wrote that correctly?
You don't just have a singular point, you have a singular plane $x_1 = x_2$. When you start close to a point on this plane (with $x_2 \ne 0$), $\dot{x_1}$ will be very large, either positive or negative depending on the sign of $x_2/(x_2 - x_1)$. Unfortunately, if $x_2 > 0$ the sign of this derivative is opposite to the sign of $x_1 - x_2$, thus $|x_1 - x_2|$ will decrease: solutions will head into the singular plane and cease to exist. There are no solutions that leave the neighbourhood of $(1,1,1)$.