# 2nd-order backward difference approximation for temporal terms of Euler–Bernoulli beam

I am trying to solve the Euler–Bernoulli beam numerically in structural dynamics analysis: $$\frac{\partial^2w(x,t)}{\partial t^2}= \frac {q(x,t)}{\mu(x)}-\frac {1}{\mu(x)}\frac {\partial^2}{\partial x^2}\left[E(x)I(x)\frac {\partial^2w(x,t)}{\partial x^2}\right]-\beta(x) \frac {\partial w(x,t)}{\partial t}$$ I use backward difference approximation for the temporal terms and central difference approximation for spatial terms. The backward difference approximations for the temporal terms referenced to this: $$\frac{\partial^2w(x_i,t_{j+1})}{\partial t^2} = \frac{2w(x_i,t_{j+1})-5w(x_i,t_j)+4w(x_i,t_{j-1})-w(x_i,t_{j-2})}{\Delta t^3}$$ $$\frac{\partial w(x_i,t_{j+1})}{\partial t} = \frac{3w(x_i,t_{j+1})-4w(x_i,t_j)+w(x_i,t_{j-1})}{2\Delta t}$$ The discretion for spatial terms is skipped here, and the term $$\frac {q(x,t)}{\mu(x)}$$ will be $$\frac {q(x_i,t_{j+1})}{\mu(x)}$$. Sort out those terms and organize them as $$Ax=B$$ and I simply write them in this form: $$[A][w(t_{j+1})]=[B(t_j,t_{j-1},t_{j-2},q(t_{j+1}))]$$ where A is the matrix of coefficients, $$w(t_{j+1})$$ is the vector of $$w$$ at $$j+1$$ time step along the beam, and $$B$$ is the vector of the combination of some known values and unknown values.

My question is: the RHS of this system has $$q(x_i,t_{j+1})$$ which are values that come from the next time step, and it is not known at the current time step, therefore, how to deal with these unknown values when solving this equation of system? I am wondering whether the predictor-corrector method can deal with this problem, or if another method will be more suitable.

Additional information: I am doing structural dynamics simulation, so the load is calculated from another part of the code every time step, so I cannot know the next time step's load in advance. And using backward differences for temporal terms because of the stability consideration.

• This is an implicit method, meaning you have to solve a linear system to obtain the values at the current time step. To start the iteration, use an explicit method with high order of accuracy. Jan 26 at 18:12
• @whpowell96, thanks a lot, but I am still confused, what's the meaning of "To start the iteration, use an explicit method with high order of accuracy.", if I start the iteration with the explicit method, how to shift to an implicit method? Jan 26 at 18:40
• Use an explicit method to compute the first three solution values then use this method to compute all later values. Jan 26 at 20:29
• @whpowell96, thanks, it seems that your comment is to address the initializing of the implicit method. But my problem is how to deal with the unknown values 'q(j+1)' of the next time step in RHS of the linear system. Jan 26 at 21:13
• How is $q$ not known here? Usually $q$ refers to prescribed load. If $q$ depends on $w$ somehow so it is not known, then you have to solve an implicit equation for $w(t_{j+1})$ Jan 26 at 21:34