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.

Improve this question
$\endgroup$
5
  • $\begingroup$ I think that it would be simpler if you first rewrite your equation in nondimensional form. $\endgroup$
    – nicoguaro ♦
    Commented Aug 20, 2018 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$
    – GGG
    Commented Aug 20, 2018 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
    Commented Aug 21, 2018 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$
    – GGG
    Commented Aug 21, 2018 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$
    – GGG
    Commented Aug 21, 2018 at 15:07

1 Answer 1

Reset to default
2
$\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.

Improve this answer
$\endgroup$
7
  • $\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$
    – GGG
    Commented Aug 22, 2018 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
    Commented Aug 23, 2018 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
    Commented Aug 23, 2018 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
    Commented Aug 23, 2018 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
    Commented Aug 23, 2018 at 18:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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