# Coupled differential equations into system of first-order equations implicitly

I am looking to solve the following equations numerically:

$$a x=\frac{d}{dt}\left(f(x,y,t)\frac{dy}{dt}\right),\quad b y=\frac{d}{dt}\left(g(x,y,t)\frac{dx}{dt}\right)$$

For arbitrary functions $$f$$ and $$g$$ and constants $$a$$ and $$b$$. I am struggling to find a way to transform this into a system of first order differential equations that I can pass into a solver. It looks like I will need to define these implicitly, but I'm not sure how to do that.

My best attempt so far is the following:

\begin{align}z_1&=f(x,y,t)\frac{dy}{dt}\\ z_2&=g(x,y,t)\frac{dx}{dt}\\ z_3&=ax\\ z_4&=by\\ \end{align}

$$$$\begin{pmatrix}z_1\\z_2\\z_3\\z_4\end{pmatrix}'=\begin{pmatrix}0&0&1&0\\0&0&0&1\\0&\frac{a}{g(t,x,y)}&0&0\\\frac{b}{f(t,x,y)}&0&0&0\end{pmatrix}\begin{pmatrix}z_1\\z_2\\z_3\\z_4\end{pmatrix}$$$$

However, this seems fairly inelegant and assumes that you are always able to divide by $$f$$ and $$g$$. I'm trying to keep this as general as possible, so don't want to make that assumption. Is there a better way to turn this into a system of first order differential equations implicitly? Thanks!

• You always have to exclude $f=0$ and $g=0$ from the domain of the ODE system, as on those surfaces the order of the system collapses, the system becomes singular there. Now you could ask under what circumstances a solution can be prolonged into these singular sets and continued (uniquely?) on the other side. – LutzL Feb 21 at 9:57

I would rather define $$u = x'$$ and $$v = y'$$. Then, your system becomes $$\begin{equation*} \begin{pmatrix} x'\\y'\\fv'\\gu' \end{pmatrix} = \begin{pmatrix} u\\v\\ ax -f_xuv -f_yv^2\\ by -g_xu^2 - g_yuv \end{pmatrix} \end{equation*}$$
As long as $$f\neq 0\neq g$$, this is a nonlinear system of first order ODEs. If, at a certain point, $$f=0$$, or $$g=0$$, or both, then this system is a differential-algebraic equation (DAE, see here and here). Numerical methods for DAEs are well studied and available. For instance, see Matlab built-in function.