1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I think that it would be simpler if you first rewrite your equation in nondimensional form. $\endgroup$ – nicoguaro Aug 20 '18 at 22:20
  • $\begingroup$ @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. $\endgroup$ – Mehrdad Yousefi Aug 20 '18 at 22:38
  • 1
    $\begingroup$ 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? $\endgroup$ – David Ketcheson Aug 21 '18 at 6:08
  • $\begingroup$ @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. $\endgroup$ – Mehrdad Yousefi Aug 21 '18 at 14:46
  • $\begingroup$ @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$. $\endgroup$ – Mehrdad Yousefi Aug 21 '18 at 15:07
1
$\begingroup$

Backward Euler and the implicit trapezoidal rule are both unconditionally stable for this problem. If you're seeing instability then you haven't implemented them correctly.

$\endgroup$
  • $\begingroup$ I checked my implementation by comparing to analytical solution and that's fine. But, the problem is because this ODE is a part of bigger simulation, I reckon it become unstable because of my Navier-Stokes solver not this ODE. Is it possible to show their unconditional stability mathematically? I mean by using something like Von-Neumann stability analysis perhaps. $\endgroup$ – Mehrdad Yousefi Aug 22 '18 at 12:55
  • $\begingroup$ Even if they are unconditionally stable, in an actual implementation you don't have true unconditional stability since it can be limited by the convergence of Newton's method. Production codes always have to do globalization techniques to the nonlinear solve, along with sophisticated extrapolators for initial point predictions, in order to be truly unconditionally stable. Even then, adaptivity of time stepping is necessary to pull back time steps when Newton's method diverges. There's a ton of tricks and heuristics to make this work, which is why I'd recommend a production ODE solver call. $\endgroup$ – Chris Rackauckas Aug 23 '18 at 13:50
  • $\begingroup$ @ChrisRackauckas for this equation you don't need Newton. Since it's a scalar, linear problem even implicit methods can be explicitly. $\endgroup$ – David Ketcheson Aug 23 '18 at 18:53
  • $\begingroup$ @MehrdadYousefi It sounds like you are actually solving a problem very different from the one you described. You should open a new question with the details of what you're really trying to do. $\endgroup$ – David Ketcheson Aug 23 '18 at 18:54
  • $\begingroup$ "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. $\endgroup$ – Chris Rackauckas Aug 23 '18 at 18:56

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.