0

I have the following generalized system of differential equations:

$$ \left( \frac{du}{dt}, \frac{dh}{dt}, \frac{dv}{dt}, \frac{dy}{dt} \right) = \left( f_1(u, h, a), f_2(u, h), f_3(u, h, u'), v \right) $$

where $a = \frac{dv}{dt}$, $u' = \frac{du}{dt}$, and $h' = \frac{dh}{dt}$

Given $ (u, h, v, y) = (0, h_0, v_0, y_0) $ and $ (u', h', a, v) = (0, 0, a_0, v_0) $ at $t = 0$.

All DEs except $\frac{dy}{dt}$ are non-linear but autonomous since they have no explicit time dependence.

How does one go about solving this system?

My attempts have been to use the Forward Euler method to solve the equations simultaneously which gives some physical solution. However my attempts to introduce the classic Runge-Kutta 4th order method was unsuccessful -- the solution was unphysical and unstable. I suspect this due to the derivative dependence of the $f_1$ and $f_3$ functions. My next idea was to simulataneously solve these two functions and eliminate this dependence on the derivatives. I then plugged these reduced functions into the solver, which gave rise to a stable solution however the second boundary condition $ (u, h, v, y) = (0, 0, v_{max}, y_{max}) $, I am not sure if the decoupling of the differential equation has been done conservatively i.e. if the original system and the reduced system still give the same solution.

I am not deeply familiar with the depths of numerical analysis however I am quite confident that I am not using the right tool for this job, am I missing something here?

Edit: I am attaching the forms of these equations below.

$f_1$ was derived based on the Transient Bernoulli's equation, $f_2$ was derived due to a system of constraints concerning the flow of fluid, $f_3$ derived based on the Newton's Second Law for variable mass systems, and $f_4$'s forms follows from definition of velocity.

$$ B(h) \frac{du}{dt} + C(h) {u}^2 + D(h) \frac{P_1}{\rho_w} - \frac{P_{atm}}{\rho_w} + (\frac{dv}{dt} + g)h = 0 $$

$$ \frac{dh}{dt} = \frac{A_{out}u}{A(h)} $$

$$ m(h) \frac{dv}{dt} = \rho_w A_{out} u^2 - \rho_w A \left( h \frac{du}{dt} + \frac{A}{A(h)} u^2 \right) - C_d v - m(h)g $$

$$ \frac{dy}{dt} = v $$

Where, $A(h), B(h), C(h)$ and $D(h)$ are functions of h which do not have a strict explicit forms ($B(h)$ is an integral and $A(h), B(h), C(h)$ are step-functions). All the other symbols are constants.

Source: https://www.et.byu.edu/~wheeler/benchtop/pix/thrust_eqns.pdf, the Water Impulse Phase

Prajval K
  • 89
  • 5
  • You already reported the obvious. Now it is time to go one step back. It seems that there was previous a higher-order system that you transformed into a first-order system. Please tell the structure of this original system and how you motivate the steps towards the first-order system. – Lutz Lehmann May 04 '23 at 09:08
  • @LutzLehmann I am not entirely sure of what you are asking, but I have attached the forms of the functions and it's origins in my edited post. – Prajval K May 05 '23 at 13:53
  • Can you provide more detailed comments on your numerical tests? There are many ways in which numerically solving ODEs can become unstable and some are easily fixable. – whpowell96 May 05 '23 at 14:12
  • 1
    At a glance, it seems odd that you have both $u'$ and $v'$ in both equations 2 and 4. You should try to substitute and rearrange these equations so that $u'$ and $v'$ can be written as functions of $(u,h,v,y)$ – whpowell96 May 05 '23 at 14:16
  • @whpowell96 I have applied RK4 classic method on solving these algorithms numerically, by unstable I mean it results in u-value (which represents the rate of propellant expulsion) in the order of 10 to the 11th powers, close to 10^4 powers of velocity. In the second comment, yeah you are absolutely right, this choice of representation by the author seems odd, and as mentioned in my post, I have attempted to substitute for $u'$ and $v'$ and got equations in the form of $(u, h, v, y)$, these equations are unstable under RK4 and SSPRK3 for time step values less than 10^-4. – Prajval K May 05 '23 at 15:50
  • I would like to know if there are methods possible that are not this expensive (the current performing algorithm SSPRK3 takes 2min for this phase, since I am also using MPFR C++ bindings to reduce numerical error while solving for at least 16 digits of floating point accuracy). Also, as mentioned in the post, is there a method for solving the equations without too much manipulation so I can no if the solution of the reduced form is different from the original system given by the author. – Prajval K May 05 '23 at 15:54
  • Instead of reinventing the wheel I suggest you use some ODE solvers already implemented in libraries in your language of choice. These have many features and failsafes that eliminate many of these headaches – whpowell96 May 05 '23 at 17:58

0 Answers0