I am trying to solve $M\ddot{u}=-Ku+F_\text{ext}$ for a 2D linear elastic model with $M$ be the mass matrix,$K$ the stiffness matrix and $F_\text{ext}$ the external load vector coming from a uniformly distributed load acting on one edge of the model.(Note: $F_\text{ext}$ is not time-dependent). An explicit time-scheme is used and more specific Forward-Euler scheme. The steps that I follow are:

  1. Initial conditions $\dot{u}_0=0$ $u_0=0$
  2. Solve $M\ddot{u}_n=-Ku_n+F_\text{ext}$ using an iterative solver
  3. Update $u_\text{n+1}=u_n+dt\dot{u}_n$
  4. Update $\dot{u}_\text{n+1}=\dot{u}_n+dt\ddot{u}_n$
  5. Go back to 2 for next time step

Based on this implementation I noticed that the output values (valocity,displacment,acceleration) go to infinity.What is the main issue that can cause this problematic behaviour?I want to note that the used time-step is small $10^{-6}$ so I don't think is a stability issue. Here is the main routine:

for(int i=0;i<2*NN;i++){

for(int i=0;i<2*NN;i++){
    double sum=0;
    for(int j=0;j<2*NN;j++){


BoundaryCondForRHS(NN,NEy,dbc,RHS);//rows connected with BC are set to zero

ConjugateGradient(2*NN,M_global,RHS,accel);//find acceleration at t->n

for(int i=0;i<2*NN;i++){
   displ[i]=dt*veloc[i]+displ[i]; //displ at t->n+1
   veloc[i]=dt*accel[i]+veloc[i]; //veloc at t->n+1
  • $\begingroup$ Do you mean Backward Euler? Forward Euler is not implicit. $\endgroup$ Oct 12, 2020 at 22:02
  • $\begingroup$ How are you coupling $\ddot u$ and $u$? It seems like you are actually solving $M \ddot u_{n+1}=K u_n+ F_{ext}$ $\endgroup$
    – Charlie S
    Oct 12, 2020 at 22:11
  • $\begingroup$ @StevenRoberts I meant explicit toy are right $\endgroup$
    – spyros
    Oct 12, 2020 at 22:12
  • $\begingroup$ @CharlieS I am using $\ddot{u}_n$ to update velocities.Is my implementation wrong? $\endgroup$
    – spyros
    Oct 12, 2020 at 22:14
  • $\begingroup$ Why not use backward Euler? You'll have to do a linear solve either way. Might as well get better stability. $\endgroup$ Oct 12, 2020 at 22:17

1 Answer 1


Your step two is to solve the original ODE, which doesn't make sense. I'll write out the steps for applying Forward Euler to your second order ODE. Forward Euler solves the first order ODE $$ M \dot{y} = f(y) $$ with the steps $$ \begin{align} M k_1 &= f(y_n) \\ y_{n+1} &= y_n + dt \, k_1 \end{align} $$

Let $v = \dot{u}$. You finite element problem can be rewritten as $$ \begin{bmatrix} I & 0 \\ 0 & M \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}' = \begin{bmatrix} v \\ K u + F_{ext} \end{bmatrix}. $$

When we apply backward Euler to this we get $$ \begin{align} \begin{bmatrix} k_1 \\ M \ell_1 \end{bmatrix} &= \begin{bmatrix} v_n \\ K u_n + F_{ext} \end{bmatrix}, \\ \begin{bmatrix} u_{n+1} \\ v_{n+1} \end{bmatrix} &= \begin{bmatrix} u_n + dt \, k_1 \\ v_n + dt \, \ell_1 \end{bmatrix} \end{align} $$

So the steps are

  1. Initial conditions $u_0$ and $\dot{u}_0$
  2. Solve the linear system $M \ell_1 = K u_n + F_{ext}$ for $\ell_1$
  3. $u_{n+1} = u_n + dt \, \dot{u}_n$
  4. $\dot{u}_{n+1} = \dot{u}_n + dt \, \ell_1$
  5. Go back to two

Edit: Now that the code is posted, it looks like you are following these steps. I suspect if you look at the eigenvalues of the Jacobain $$ \begin{bmatrix} 0 & I \\ M^{-1} K & 0 \end{bmatrix}, $$ there will be some (if not all) along the imaginary axis, which is outside the stability region. I would recommend trying Backward Euler or the Newmark method.

  • $\begingroup$ I think you wrote exactly what I did ,with $l1 = \ddot{u}_n$ $\endgroup$
    – spyros
    Oct 12, 2020 at 22:34
  • $\begingroup$ So to be crystal clear ,the problem lies to the time-integration method and not to implementation? $\endgroup$
    – spyros
    Oct 12, 2020 at 22:49
  • $\begingroup$ Yes, it looks like the steps are the same. I got a bit confused with the $\ddot{u}_n$ notation. I don't see any issues with the code. For a short timespan, I would expect forward Euler to converge, even with eigenvalues outside the stability region. If possible, I would try to construct the matrix above and check the eigenvalues. $\endgroup$ Oct 12, 2020 at 23:08
  • $\begingroup$ Is it possible, that the time step is not small enough?For example may a time step of $10^{-9}$ or something like that is needed $\endgroup$
    – spyros
    Oct 12, 2020 at 23:33
  • 1
    $\begingroup$ Your linear elastic problem will oscillate and thus have imaginary eigenvalues. As stated above, for purely imaginary eigenvalues, explicit Euler is unconditionally unstable: that is, no matter how small your time step, your solution will diverge. If you make the time step small, it will grow slower but it will still grow. You really need to use a different integration method like Newmark, trapezoidal rule or something which is stable for imaginary eigenvalues. $\endgroup$
    – Daniel
    Oct 13, 2020 at 7:05

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.