Second Order Finite Difference Scheme to Solve Initial Value Problem

Question

I would appreciate help in a problem that I cannot figure out where am I doing the mistake. The mistake could be something I am unaware of so please excuse my lack of knowledge.

I am trying to solve the equation of motion, $$M \ddot{x} + K x = F,$$ where $M$, $K$, and $F$ are variables independent of $x$, and $\ddot{x}=\frac{d^2x}{dt^2}$. Using a central difference finite scheme, I rewrote the equation as,

$$M \frac{ x_{i+1} - 2x_i + x_{i-1} }{ {\Delta t} ^2 } + Kx = F.$$

Then, I am solving this equation for $x_{i+1}$ at every time step $i$ as,

$$ x_{i+1} = \frac{{\Delta t}^2}{M} \cdot \left( F - K x + \frac{M}{{\Delta t}^2} \cdot 2 x_i - \frac{M}{{\Delta t}^2} \cdot x_{i-1} \right).$$

Now to solve this problem, one needs two initial conditions, to my knowledge, namely, $x(0)$ and $\dot{x}(0)$. At the first time step, I am setting $x_i$ to the $x(0)$ initial condition. Now I am stuck with two problems:

• What value should I set $x_{i-1}$ at first time step?
• How can I apply the initial condition $\dot{x}(0)$?

Answer 1

One trick is to introduce "ghost cells" like $x_{-1}$ with the definition per boundary condition $$ \frac{x_1-x_{-1}}{2Δt}=\dot x(0)\implies x_{-1}=x_1-2Δt\,\dot x(0) $$ Now you can either add this to the system, or directly eliminate $x_{-1}$ against the differential equation centered at $x_0$.

Another way to use the boundary condition with a sufficiently high order is to use a forward differentiation formula of that order, like in the second order $$ \frac{-3x_0+4x_1-x_2}{2Δt}=\dot x(0). $$

Answer 2

$x_{i-1}=x_{i}-\dot x(0)\Delta t$

Answer 3

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

• In the firt step,lets move everything to the adimensional form: \begin{align} & x = {F \over k}\,X,\quad t = \root{M \over k}T \\[2mm] & \substack{\mbox{which yields, with} \\[1mm] \ds{\dot{a} \equiv \totald{a}{t},\quad b' = \totald{b}{T}}} \quad \left\{\begin{array}{l} \ds{X'' + X = 1} \\[2mm] \ds{X\pars{0} = {k \over F}\,x\pars{0} \equiv X_{0}} \\[2mm] \ds{X'\pars{0} = {\root{kM} \over F}\,\dot{x}\pars{0} \equiv V_{0}} \end{array}\right. \end{align}
• With $\ds{Y = X'}$, I'll have the couple of equations \begin{align} &\left.\begin{array}{rcl} \ds{X'} & \ds{=} & \ds{Y} \\ \ds{Y'} & \ds{=} & \ds{1 - X} \end{array}\right\}\quad\mbox{with}\quad \left\{\begin{array}{rcl} \ds{X\pars{0}} & \ds{=} & \ds{X_{0}} \\ \ds{Y\pars{0}} & \ds{=} & \ds{V_{0}} \end{array}\right. \end{align} and the "discrete version" $$ \left\{\begin{array}{rcl} \ds{X_{n + 1}} & \ds{=} & \ds{X_{n} + Y_{n}\,\Delta T} \\ \ds{Y_{n + 1}} & \ds{=} & \ds{Y_{n} + \pars{1 - X_{n}}\Delta T} \end{array}\right. $$

---

With $\ds{\color{red}{X_{0}} = \color{blue}{V_{0}} = 0}$, the following picture illustrates $\ds{10000}$ iterations of $\ds{\color{red}{X\pars{T}}}$ and $\ds{\color{blue}{Y\pars{T} = X'\pars{T}}}$ with $\ds{\Delta T = 0.005}$: