Polar coordinates for vector field to find sticking flow

Question

I am currently working on an impacting system which is basically just a spring damper and a circular enclosure. Because of the rotational symmetry of the problem I need the vector field in polar coordinates to derive the sticking flow along the enclosure. It should be an easy problem, but contrary to my expectation there are no scores of textbooks which contain the solution and I am not used to mechanical problems.

The system in Cartesian coordinates is the following: $$ \left( \begin{array}{c} \dot{x} \\ \ddot{x} \\ \dot{y} \\ \ddot{y} \\ \end{array} \right)= \left( \begin{array}{cccc} 0 & 1 & 0 & 0 \\ -\nu^2 & -\gamma & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -\nu^2 & -\gamma \\ \end{array} \right) \left( \begin{array}{c} x \\ \dot{x} \\ y \\ \dot{y} \\ \end{array} \right) $$

I have replaced $x$, $y$ and their derivatives by $x=r \cos(\theta)$, $y=r \sin(\theta)$ and their derivatives as was also done in this question. Solving for $\ddot{r}$ and $\ddot{\theta}$ I obtained $$ \ddot{r}= -r \nu ^2-\gamma \dot{r}+r \dot{\theta} ^2 $$ and $$ \ddot{\theta}=-\frac{\left(r \gamma +2 \dot{r}\right) \dot{\theta}}{r} $$ which I verified with mathematica.

But when I test using the following Matlab-code

nu2=1.96; gam=0.28;
ODE1=@(t,z) [0 1 0 0;-nu2 -gam 0 0;0 0 0 1;0 0 -nu2 -gam]*z;
ODE2=@(t,z) [z(2); 
            -nu2*z(1)-gam*z(2)+z(1)*z(4)^2; 
             z(4); 
            -z(4)*(gam*z(1)+2*z(2))/z(1)];
z0=[1;0;1;0.5];
tspan=[0, 1];
z0p=[1.4142 0.3536 0.7854 0.3536];
[T1,Z1]=ode45(ODE1,tspan,z0);
[T2,Z2]=ode45(ODE2,tspan,z0p);
Z2c=[Z2(:,1).*cos(Z2(:,3)),Z2(:,1).*sin(Z2(:,3))];
hold on
plot(Z1(:,1),Z1(:,3));
plot(Z2c(:,1),Z2c(:,2));
hold off
legend('Cartesian','Polar')

I find that they are not equivalent. Now I am wondering whether there is a conceptual error as I am fairly confident in my derivation. So I am starting to think that my approach must be wrong and I would happy if someone could tell me where my error lies.

Update

Found the problem, answer below.

Answer 1

It turns out the problem was in the calculation of the initial velocities which should be $$ \dot{r}=\dot{x}\cos(\theta)+\dot{y}\cos(\theta) \quad \text{and} \quad\dot{\theta}=\frac{\dot{y}\cos(\theta)-\dot{x}\cos(\theta)}{r} $$ Then the the system above is correct and both formulations yields the same plot of the correct initial conditions in polar coordinates are given which are

z0p=[1.4142    0.3536    0.7854    0.2500]