I'm trying to solve the fourth order ODE characterising elastic curves to compute some elastic curves:
$$ \gamma'''' + 3 \langle \gamma'', \gamma''' \rangle \gamma' + \frac{3}{2} \langle \gamma'',\gamma'' \rangle \gamma'' = 0$$
First of all I reformulated the fourth order ODE as first order ODE:
$$1)\ a_1 = \gamma'$$
$$2)\ a_2 = \gamma'' = a_1'$$
$$3)\ a_3 = \gamma''' = a_2'$$
$$4)\ a_3' = \gamma'''' = -3 \langle a_2, a_3 \rangle a_1 - \frac{3}{2} \langle a_2,a_2 \rangle a_2$$
I want to compute planar curves, i.e. $\gamma(t) = (\gamma_x(t),\gamma_y(t))$ for $t \in [a,b]$. I get this "second" system of ODEs, namely one for each dimension, since the $a_i = (a_{i,x},a_{i,y})$ are two dim vector maps:
$$1i)\ a_{1,x} = \gamma'_x$$
$$1ii)\ a_{1,y} = \gamma'_y$$
$$2i)\ a_{2,x} = \gamma''_x$$
$$2ii)\ a_{2,y} = \gamma''_y$$
$$3i)\ a_{3,x} = \gamma'''_x $$
$$3ii)\ a_{3,y} = \gamma'''_y$$
$$4i)\ a_{3,x}' = \gamma''''_x = -3 \langle a_2, a_3 \rangle a_{1,x} - \frac{3}{2} \langle a_2,a_2 \rangle a_{2,x}$$
$$4ii)\ a_{3,y}' = \gamma''''_y = -3 \langle a_2, a_3 \rangle a_{1,y} - \frac{3}{2} \langle a_2,a_2 \rangle a_{2,y}$$
However when I implement this in python and let it solve(e.g I tried scipy.integrate.odeint and scipy.integrate.solve_ivp) the results don't look very promising. I played around with the eight initial values without improvement.
I think the problem are the last two equations. I'm not very familiar with solving ODEs, but in 7) and 8) I have not only the component-wise functions but because of the scalar product both, x and y components:
$$4ii)\ a_{3,y}' = \gamma''''_y = -3( a_{2,x}a_{3,x}+a_{2,y}a_{3,y}) a_{1,y} - \frac{3}{2} (a_{2,x}a_{2,x}+a_{2,y}a_{2,y}) a_{2,y}$$
Do the ODE solver not like this "cross reference" of the functions? But I don't know why this should be a problem and what to change...
Is this even the right way to tackle this problem?
Please let me know if something is not clear and thanks for anyone who can put me in the right direction!