First published in internet mathematics forums

Posted: Friday 4th August 2017

I derived this system of differential equations this week as I was researching possible profiles for water diffusers & nozzles, used when joining pipes with different bores.

**Mathematics:*** Solve this system of differential equations.*

\[ x' = y^{-2} \]

\[ y' = - \sqrt {(t+1)^2- y^{-4}} \]

\(x'\) and \(y'\) are derivatives with respect to \(t\).

I have obtained a **numerical** solution (which was non-trivial because of the numerical instability of the Euler method with this system of differential equations) but I am curious to know** "does an analytical solution exist?"**, which would be more efficient and convenient to use.

**Derivation of the system of differential equations**

Water is accelerated in a nozzle or a pipe of reducing width, which is rotationally symmetrical about the X-axis, with the bore, the inner diameter and the inner radius proportional to a function \(y(x)\) .

Neglecting viscosity and considering averages for simplicity, the velocity of the water in the X-axis is inversely proportional to the bore's area of cross section and to \(y^2\).

\[x' = y^{-2}\]

The average velocity \(r'\), of a radial element, a thin slice of water "pie", is composed of the vector addition of the velocity along the X-axis \(x'\) and in the radial direction \(y'\), which are related in magnitude by Pythagoras,

\[ r'^2 = x'^2 + y'^2 \]

So

\[ y' = \sqrt{r'^2 - x'^2} \]

Assume an acceleration \( r' = at + u \), at a time \(t\), with acceleration a and initial velocity u, but for simplicity here, both a and u are assumed to be 1.

\[ r' = t + 1 \]

Therefore

\[ y' = - \sqrt{ (t+1)^2 - y^{-4}} \]

choosing the negative root corresponding to an radially inward \(y'\) when water accelerates in a nozzle.

**Numerical Solution**

The numerical instability was managed by writing a computer program which could calculate in t-increments corresponding to the square root of linear increments in \(t^2\).

\( t_{i+1} = t_i + \Delta t_i \) where \( \Delta t_i = \min(h_1,\sqrt{t_i^2+h_2}-t_i) \) and \(h_1 \), \( h_2 \) are step size constants.

**Approximate solution for \( (x,y,t) = (0,1,0)\)**

With the initial conditions \( t=0, x=0, y=1 \) then \( x' = r' = 1 \) and \( y'=0 \).

Assuming that then \( y' << x' \) for all \( t>0 \) then approximately

\[ x' = t + 1 \]

Integrating with respect to t and substituting for t (or simplifying this equation for linear acceleration, \( x'^2 - u^2 = 2ax \) with \(a=1\) and \(u=1\) ) gives

\[ x'^2 - 1 = 2x \]

Substituting for \( x' = \sqrt{2x+1} \) in the system equation \( x' = y^{-2} \) and rearranging for y gives

\[ y = (2x+1)^{-0.25} \]

As this graph shows, this is a good approximate solution for these starting conditions.

Note on numerical instability

The graph also shows how numerical stability interrupted the numerical solution part plotted using the free on-line Two Dimensional Differential Equation Solver and Grapher V 1.0. Investigate the numerical instability by selecting the option *"System of first order DEs: x' = f(x, y, t), y' = g(x, y, t)" *and typing for x', y' -

x' = y^(-2)

y' = -1*sqrt((t+1)^2-y^(-4))

with initial values (x=0, y \(\ge\)1)