what is the limit of $\theta$ if $r\to 0$ when we have these conditions? the problem is based on this picture. 
at beginning or we say $t=0$, $P$ is a circle of which the center is at the point $(0,r)$, $r_0=1$ is the initial radius of this circle. $AB$ is a vector which has an angle $\theta$ from the $x$ axis, and this vector $AB$ will always on the tangent of the circle (the side close to $O$). We suppose that $A_0=(x_0,y_0)$ and the initial value of $\theta$ is $\theta_0$. we have also $|AP|=l=2$. the angle betwen $PA$ and $x$ axis is noted $\varphi$.
The vector $AB$ moves on the tangent of the circle with a velocity $-x$, so $A$ will be closer and closer to O. when the distance $AP=l<2r$, the radius becomes smaller, and then we have new $P,r$, but the new $r$ should satisfy $AP=2r$. 
The question is : when $A\to O,r\to 0$, what is the limit of $\theta$?
I tried to use the differential method to get the result, but didn't make it. I know that the angle $\angle PAQ$ is always $\pi/6$ and $\dot{l}=2\dot{r}$. I don't know how to continue. thanks very much for your help.
 A: Some analysis
I choose $x$ as a free parameter, and $y(x)$ as a function to describe the movement. You want $l=2r$ which you can rewrite as
\begin{align*}
x^2+(y-r)^2&=(2r)^2\\
x^2+y^2-2yr+r^2&=4r^2\\
3r^2+2yr-(x^2+y^2)&=0\\
r&=\frac{-y+\sqrt{3x^2+4y^2}}{3}
\end{align*}
The other solution of the quadratic equation would result in $r<0$ and can therefore be omitted. Now you know that
$$\overrightarrow{PA}=\begin{pmatrix}x\\y-r\end{pmatrix}
=\frac13\begin{pmatrix}3x\\4y-\sqrt{3x^2+4y^2}\end{pmatrix}$$
Rotate that right by $60°$ and scale it by $6$ and you get
$$12\overrightarrow{PQ}=
\begin{pmatrix}
1&\sqrt3\\-\sqrt3&1
\end{pmatrix}\cdot3\overrightarrow{PA}=
\begin{pmatrix}
3x+\sqrt3\left(4y-\sqrt{3x^2 + 4y^2}\right)\\
-3\sqrt{3} x + 4y - \sqrt{3x^2 + 4y^2}
\end{pmatrix}
$$
That vector has to be orthogonal to the direction of the movement, namely
$$\begin{pmatrix}1\\y'\end{pmatrix}$$
so their dot product must be zero. This will give you a relation between $y(x)$ and $y'(x)$:
$$
y'=\frac{\mathrm dy}{\mathrm dx}=
-\frac{3x+\sqrt3\left(4y-\sqrt{3x^2 + 4y^2}\right)}
{-3\sqrt{3} x + 4y - \sqrt{3x^2 + 4y^2}}
$$
Interpreting this formula
Now I don't have enough experience to solve this beast analytically (if that is at all possible), and for a numeric solution we would need precise starting conditions. But a look at the slope field might be instructive:

That does look reasonable to me.
Convegence claim
One thing worth noting is that the whole problem is scale invariant (at least if you are in the regime where the circle shrinks with the approaching point). If you multiply all lengths with a common factor, the slope remains the same. So if the slope converges towards the end, it means that the final direction must be a fixpoint, so starting from that direction at any distance will cause movement in a straight line towards the origin. But there is only one direction which matches movement directly to the origin in the initial situation, and that is horizontal. That's because the tangent at the origin is horizontal.
So I'd say the slope converges towards horizontal eventually, no matter the initial conditions. In the limit, $\theta=0$.
Numerical experiments
The convergence is very slow towards the end, though. In the following animation, I'm using the distance $OA$ as a scale factor for the speed, and dowards the end the change in the slope of the line $AB$ is hardly perceptible any more. So I can see how numeric experiments would indicate a non-zero limit. Nevertheless, the apparent final slope cannot be the real final slope, because if it were, then slopes along that line would have to point directly towards the origin, which they do not.
Also note that the convergence only appears slow in this parametrization with adaptice speed. If you use fixed speed, you will reach the point where $x=0$ in finite time. But obtaining precise numeric results for this case is hard, which is the reason why I consdier this adaptiveness and the resulting slow convergence to be reasonable.

Change in angle
In a comment, you asked about characterizing the convergence in terms of a Lyapunov exponent. For that we'd need to fix a reasonable parameter space. Since the problem is scale invariant, I'd use the angle which $OA$ makes againt horizontal as the sole description of phase space. Let's call that angle $\alpha$, with $\tan\alpha=y/x$. For the time dependence, I'll consider a point at unit distance from the origon, traveling at unit speed towards $Q$. This formulation models the rescaling I used above. Now consider the following differential equation:
\begin{align*}
\frac{\cos(\alpha+\mathrm d\alpha)}{\sin(\alpha+\mathrm d\alpha)} =
\frac{\cos\alpha-\mathrm d\alpha\sin\alpha}{\sin\alpha+\mathrm d\alpha\cos\alpha}
= \frac{x-y\,\mathrm d\alpha}{y+x\,\mathrm d\alpha}
&= \frac{x+\mathrm dt\,\Delta x}{y+\mathrm dt\,\Delta y}
\\
(x-y\,\mathrm d\alpha)(y+\mathrm dt\,\Delta y) &=
(y+x\,\mathrm d\alpha)(x+\mathrm dt\,\Delta x)
\\
x\,\Delta y\,\mathrm dt-y^2\,\mathrm d\alpha &=
y\,\Delta x\,\mathrm dt+x^2\,\mathrm d\alpha \\
(x^2+y^2)\mathrm d\alpha = \mathrm d\alpha &= (x\,\Delta y-y\,\Delta x)\mathrm dt
\\
\dot\alpha = \frac{\mathrm d\alpha}{\mathrm dt} &= x\,\Delta y-y\,\Delta x
\end{align*}
Here I assume $x=\cos\alpha,y=\sin\alpha$ and also assume $(\Delta x,\Delta y)$ to be a unit length vector pointing from $A$ to $Q$:
$$\begin{pmatrix}\Delta x\\\Delta y\end{pmatrix}=
\frac{-1}{\sqrt{1+y'^2}}\begin{pmatrix}1\\y'\end{pmatrix}$$
at least for the case where $y'>0$. For $y'<0$ we want the opposite sign; those are the parts where we move down and right. But for the simpler case of $y'<0$ we get
\begin{align*}
\dot\alpha=
\frac{\mathrm d\alpha}{\mathrm dt}&=
\frac{y-xy'}{\sqrt{1+y'^2}}
\end{align*}
Fixing the sign over the whole range, you get the following as a plot of the change in angle (or rather its negative, for better placement of labels):

Or here a close-up of the last degree:

As you can see, the smaller the angle, the smaller its change will be. But the change will always be negative, so the tangent will always move towards horizontal.
Now that I think about it, it will be difficult to express this in terms of a (global) Lyapunov exponent: On the one hand there is a high change rate for large angles, whereas for arbitrarily small angles, the change rate can become arbitrarily slow as well. But I've not much experience dealing with Lyapunov exponents, so perhaps you can find a better formulation for this.
Looking at these graphs, it should be possible to show that $\dot\alpha<0$ for $\alpha\in(0°,90°]$, thereby demonstrating the convergence
$$\lim_{t\to\infty}\alpha(t)=0$$
which in turn implies $\theta=0$ and $\varphi=-\frac\pi6=-30°$ in that limit.
Constant speed
Also please remember that the slow convergence is a result of adaptive step width. If you move at uniform speed, e.g. from $x_0=y_0=1$, then angle as a function of time will look like this:

You can obtain such a graph with reasonable precision if you do the adaptive step width approach, but keep track of the curve length so far and use that as the $t$ coordinate for the data point. This will give you a lot of points on the near vertical axis at the far right.
The near-vertical slope there will tell you that the difference between correctly representing the convergence and overshooting the $x=0$ situation is very slim.
A: Starting from the beautifull result of MvG, it is possible to integrate $$
y'=\frac{\mathrm dy}{\mathrm dx}=
-\frac{3x+\sqrt3\left(4y-\sqrt{3x^2 + 4y^2}\right)}
{-3\sqrt{3} x + 4y - \sqrt{3x^2 + 4y^2}}
$$
and finally obtain an equivalent  $-2\sqrt(3)\frac{x}{\ln(x)}$ for the function $y(x)$. This proves that $y(x)$ tends to $0$ for $x$ tending to $0$

The preceeding result leads to the equation of the family of curves expressed on parametric form (the parameter is $z$).
This shows that all the curves are homothetic. The factor $C$ is determined by the initial point. So, it is sufficient to draw only one curve (for example $C=1$) to show the behaviour of all. 
For small values of the parameter $z$, the point $(x,y)$ is closed to $0$. In order to clearly see the decreassing of $\theta$ while $z$ tends to $0$, it is necessary to go to very small values. This is due to the logarithmic term in the equivalent function. This is clear below on the several graphs with different magnitudes.

When $z$ tends to infinity, $x$ tends to $0$ and $y$ tends to $y_{max} = e^{-2+\Re(\tanh^{-1}(2) ) } $ approximately $= 0.234408$ where $\Re(\tanh^{-1}(2))$ is the real part of $\tanh^{-1}(2)$ which is complex.
