# Stability Criteria for Numerical Solution of Windkessel Ordinary Differential Equation

I'm trying to solve this equation (Windkessel equation) numerically as:

$$C \frac{d P}{d t} + \frac{P}{R} = Q(t)$$

Where $C$ is compliance, $R$ is resistance, $P$ is pressure, and $Q(t)$ is a known input signal (i.e. volumetric flux). I tried explicit (i.e. Forward Euler Method) as well as semi-implicit (i.e. Backward Euler Method) and trapezoidal time integration techniques. But, all of these three well-known methods fail at some values for $C$ and $R$.

First, I wanted to know how can I map the stability region of this ODE based on $R$ and $C$ values? Also, what's your suggestion in terms of the time integration method? For example, will the 4th order Runge-Kutta method help me to reach a broader stability zone or not? Any comment or idea is appreciated.

• I think that it would be simpler if you first rewrite your equation in nondimensional form. – nicoguaro Aug 20 '18 at 22:20
• @nicoguaro Its nondimensional form is: $\frac{d \hat{P}}{d \hat{t}} + \frac{\hat{P}}{\hat{R}\hat{C}} = \frac{\hat{Q}(\hat{t})}{\hat{C}}$. I know it really depends on $\hat{R}\hat{C}$, which could be assumed as a relaxation time. Any comment or suggestion is appreciated. – Mehrdad Yousefi Aug 20 '18 at 22:38
• Are C and R constants? Are they positive or negative? What do you mean when you say that backward Euler is semi-implicit? And what do you mean when you say the methods fail? – David Ketcheson Aug 21 '18 at 6:08
• @DavidKetcheson Both $R$ and $C$ are positive constant real values. My purpose from semi-implicit (i.e. Backward Euler Method) is that, in fact, I discretized the above equation as: $C \frac{P(t+\Delta t) - P(t)}{\Delta t} + \frac{P(t+\Delta t)}{R} = Q(t)$. I call it semi-implicit because I just used the the next time step value for $\frac{P(t+\Delta t)}{R}$ but because the $Q(t+\Delta t)$ is unknown for me, I approximate it by $Q(t)$ in the current time step. My purpose from failing is that they become unstable and produce wrong results. – Mehrdad Yousefi Aug 21 '18 at 14:46
• @DavidKetcheson In fact, this equation is a small part of a bigger simulation. As a result, $Q(t+\Delta t)$ is unknown for me because I calculate $Q$ from Navier-Stokes equation and update my pressure according to this ODE. But, I'm trying to find a A-stable or L-stable time integrator which at least could guarantee unconditional stability for me at some combination of $R$ and $C$. – Mehrdad Yousefi Aug 21 '18 at 15:07

• "Since it's a scalar, linear problem even implicit methods can be explicitly.", oh true. For some reason I though $Q=Q(t,P)$ and nonlinear, my bad. – Chris Rackauckas Aug 23 '18 at 18:56