# Resolving a stiff hyperbolic problem with Neumann boundary conditions

I am trying to numerically resolve the equation for an Euler-Bernoulli beam that is inextensible, unshearable, and subject to planar deformations:

$$\rho I(s) \frac{\partial^2 \theta}{\partial t^2}(s,t) = E \frac{\partial}{\partial s} \left[I(s) \left( \frac{\partial \theta}{\partial s}(s,t) - \kappa_p(s,\theta,t) \right) \right]\,,\,\,\,\,\, t>0\,,\,\,\,\,\, 0<s<L\,,$$

$$\frac{\partial \theta}{\partial s}(0,t) = \frac{\partial \theta}{\partial s}(L,t) = 0\,,$$

$$\theta(s,0) = \frac{\partial \theta}{\partial t}(s,0) = 0\,,$$

where $\kappa_p(s,\theta,t)$ is a preferred curvature. The ends of the beam are force-free and moment-free. Here, $E \gg 1$ so that $\frac{\partial \theta}{\partial s}(s,t) \approx \kappa_p(s,\theta,t)$. I cannot directory approximate $\theta(s,t)$ with $\kappa_p(s,\theta,t)$ since $\kappa_p(s,\theta,t)$ depends on $\theta(s,t)$. I have been reading into what may be a good numerical method to use to resolve this problem, however, I am not finding very good solutions since this equation has a stiff part in both the linear term and the nonlinear forcing term. An implicit method would free me from the minisulce time-step coming from $E$, however, the phase will not be resolved properly. I can use an implicit method for the linear terms, however, using something like multi-level approximation in time for the source term also introduces errors. Also, an implicit method may not resolve the phase well if I take big time-steps. What would be the best way to resolve this problem? If the area moment of inertia, $I(s) = 0$, at the end points, would it be a good idea to simply replace the geometry of the beam where the endpoints have a finite, but small $I(s)$ so that the PDE isn't singular or is there a nice way that I can numerically treat this singularity?

Aside: To show that the system is hyperbolic, I can rewrite the PDE as a system of equations and then compute the eigenvalues of the linear part.

$$\rho I(s) \frac{\partial \theta}{\partial t} = \frac{\partial \omega}{\partial s}$$

$$\frac{\partial \omega}{\partial t} = E \left[I(s) \left( \frac{\partial \theta}{\partial s}(s,t) - \kappa_p(s,\theta,t) \right) \right]$$

Now, if $I(s)$ is strictly positive, then the eigenvalues are $\lambda_{\pm} = \pm \sqrt{E/\rho}$ which is a real quantity in this problem. Thus, the problem is hyperbolic. For non-negative $I(s)$, one recovers the same eigenvalues when taking the limit. Note that the stiffness for the problem is:

$$S = \frac{|\lambda_+|}{|\lambda_-|} = 1$$

however, depending on the numerical method, the variation in $I(s)$, and if I am not careful about the treatment of the singular points (from $I(0) = I(L) = 0$), then I will get numerical stiffness.

• Are you sure it's hyperbolic? It looks like it could also be elliptic, depending on signs. – David Ketcheson Nov 29 '15 at 18:14
• @DavidKetcheson, I should have mentioned that $\rho$, $I(s)$, and $E$ are positive. If I consider the simple case where $\rho = I(s) = 1$ and $\kappa = 0$, then I have $\theta_{tt} = E \theta_{ss}$ which is the wave equation. – namu Nov 30 '15 at 2:00
• @DavidKetcheson, it is even more clear if I write the system as $\rho I(s) \theta_t = \omega_s$ and $\omega_t = E[I(s) \theta_s - \kappa(s,\theta,t)]$ Then, the eigenvalues for the linear system are $\lambda = \pm \sqrt{E/\rho}$ and because both $E$ and $\rho$ are positive, this is a hyperbolic system. – namu Nov 30 '15 at 2:30