I am looking for a numerical method to help me solve the following system of two equations for $\theta(\xi,t)$ and $\lambda(\xi,t)$:

\begin{equation} \ddot{\theta}=\theta''+\lambda \\ \dot{\theta}^2\cos\theta+\ddot{\theta}\sin\theta=\lambda''+\lambda'\theta' \end{equation}

where $\dot{\theta} = \frac{\partial \theta}{\partial t}$ and $\theta'=\frac{\partial \theta}{\partial \xi}$.

I have the initial conditions $\theta(\xi,0)$ and $\dot{\theta}(\xi,0)$. I also have the boundary conditions $\theta(0,t)$, $\theta(\ell,t)$, $\lambda(0,t)$, and $\lambda(\ell,t)$.

I have tried to set this up for the builtin function pdepe to input into MATLAB, but I have not been able to put the equations into parabolic-elliptic form.

Thank you.

EDIT: The system stated above is a model system for me to try and get a hold of the problem. The equations I am actually trying to solve are a little different. They come from balances of linear and angular momentum for a rod which can bend. The rod has a torque, $T$, applied at its left end at time zero. I am trying to obtain the slope of the rod, $\theta$, and the internal axial and shear forces in the rod, $\lambda$ and $\mu$, for all values of space and time. The actual system is as follows:

\begin{equation} \ddot{\theta}=\theta''+\mu \\ \left( -\ddot{\theta}\sin\theta - \dot{\theta}^2\cos\theta \right)\cos \theta + \left( \ddot{\theta} \cos \theta - \dot{\theta}^2 \sin \theta \right) \sin \theta = \lambda'' - \lambda \theta'^2 - 2 \mu' \theta' -\mu \theta'' \\ \left( \ddot{\theta} \sin \theta + \dot{\theta}^2 \cos\theta \right) \sin \theta + \left( \ddot{\theta} \cos{\theta} -\dot{\theta}^2 \sin \theta \right) \cos \theta = \mu'' - \mu \theta'^2 + 2\lambda'\theta'+\lambda\theta''. \end{equation}

I have left out some of the quantities carrying physical dimension, such as the flexural rigidity and density, to focus on the mathematics of the problem.

The rod is initially straight and isn't moving. This gives me the IC's: $\theta \left( \xi,0 \right) = 0$ and $\dot{\theta} \left(\xi,0 \right)=0$.

There is a torque at the left end and no torque at the right end ($\xi = \ell$). This gives me the BC's: $\theta' \left(0^+,t\right) = T$ and $\theta' \left(\ell^-,t\right) = 0$.

The rod has no axial force at either the left or right end. This gives me the BC's: $\lambda \left( 0^+,t \right) = 0$ and $\lambda \left( \ell^-,t \right) = 0$.

The rod has no shear force at either the left or right end. This gives me the BC's: $\mu \left( 0^+,t \right) = 0$ and $\mu \left( \ell^-,t \right) = 0$.


2 Answers 2


One notes that there is no time derivative of $\lambda$. Substituting the first equation in the second one, we rewrite the system as \begin{aligned} \partial_t \theta &= \dot\theta \\ \partial_t\dot\theta &= \theta'' + \lambda \\ {\dot\theta}^2\cos\theta + (\theta'' + \lambda)\sin\theta &= \lambda'' + \lambda'\theta' \, . \end{aligned} Using finite differences on a regular grid $t_n = n\Delta t$, $\xi_i = i\Delta \xi$, one may write \begin{aligned} \theta_i^{n+1} &= \theta_i^n + \Delta t\,\dot\theta_i^n \\ \dot\theta_i^{n+1} &= \dot\theta_i^{n} + \Delta t\, \frac{\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n}{{\Delta\xi}^2} + \Delta t \, \lambda_i^n \end{aligned} and \begin{aligned} &\left(\dot\theta_i^n\right)^2\cos\theta_i^n + \left(\frac{\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n}{{\Delta\xi}^2} + \lambda_i^n\right)\sin\theta_i^n \\ =\; &\frac{\lambda_{i+1}^n - 2 \lambda_i^n + \lambda_{i-1}^n}{{\Delta\xi}^2} + \frac{\lambda_{i+1}^n - \lambda_{i-1}^n}{2\,{\Delta\xi}}\frac{\theta_{i+1}^n - \theta_{i-1}^n}{2\,{\Delta\xi}} . \end{aligned} The first two equations update $\theta$ and $\dot\theta$. The last equation provides the linear algebraic system \begin{aligned} &(1 + \tfrac{1}{4}(\theta_{i+1}^n - \theta_{i-1}^n))\lambda_{i+1}^n - (2 + {\Delta\xi}^2\sin\theta_i^n) \lambda_i^n + (1 - \tfrac{1}{4}(\theta_{i+1}^n - \theta_{i-1}^n))\lambda_{i-1}^n \\ =\; & {\Delta\xi}^2\left(\dot\theta_i^n\right)^2\cos\theta_i^n + (\theta_{i+1}^n - 2 \theta_i^n + \theta_{i-1}^n)\sin\theta_i^n \end{aligned} satisfied by $\lambda_1^n,\dots , \lambda_I^n$, which values are required in the time-update of $\dot\theta$. The initial conditions give $\theta_i^0 = \theta(\xi_i,0)$ and ${\dot\theta}_i^0 = \dot\theta(\xi_i,0)$. The boundary conditions give $\theta_0^n = \theta(0,t_n)$, $\theta_I^n = \theta(\ell,t_n)$, etc.

  • $\begingroup$ Harry, thank you so much for your reply. Please see the update to my original post for some problem background. I am going to attempt to apply finite differences to the actual problem that I am trying to solve. Let me know if you don't think this same method will work on my actual problem. $\endgroup$
    – Evan
    Commented Oct 4, 2018 at 16:22
  • $\begingroup$ @Evan it is a bit more work to solve the real problem instead of the toy-model, but a similar method may work fine $\endgroup$
    – EditPiAf
    Commented Oct 4, 2018 at 17:05

HINT: I think theres factorization you can do here...

Consider \begin{equation} \dot{\theta}^2\cos\theta+\ddot{\theta}\sin\theta \end{equation}

Then consider:

$$\frac{\partial^2 \{-\cos(f(t))\}}{\partial^2t } = f'(t)^2\cos(f(t))+f''(t)\sin(f(t))$$

Substitute $f(t)=\theta(t)$ Quite similar, no?



You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .