Consider a function $y$ with initial data $y(0)=\alpha,y'(0)=\beta$. Assume that $y$ satisfies the following equation $$y(t)=\alpha\cos(\lambda t)+\frac{\beta}{\lambda}\sin(\lambda t)-\int_0^t\frac{\sin(\lambda(t-s))}{\lambda}g(y(s))\,ds$$ where $\lambda\gg1$, $\alpha,\beta$ are given constants, and $g(y)$ is a given Lipschitz continuous function.
Choose a time step $\tau>0$ and denote $t_n=n\tau$ for $n\ge0$. Let $y^n$ be the numerical approximation of $y(t_n)$ for $n\ge0$. Based on the integral formulation, design the following time integrator via proper numerical quadratures $$ y^{n+1}=2 \cos (\lambda \tau) y^{n}-y^{n-1}-\frac{\sin (\lambda \tau)}{\lambda} g\left(y^{n}\right), \quad n \geq 1 $$ with $$ y^{0}=\alpha, \quad y^{1}=\alpha \cos (\lambda \tau)+\frac{\beta}{\lambda} \sin (\lambda \tau)-\frac{\sin (\lambda \tau)}{2 \lambda} g(\alpha) $$ Under proper stability assumption, prove the following error bound $$ \left|y\left(t_{n}\right)-y^{n}\right| \leq C \tau^{2}, \quad 0 \leq n \leq \frac{T}{\tau} $$ where $C>0$ is a constant independent of $\tau$. $$$$ From the Lipschitz condition, we see immediately that $y^1$ is enough 'close' to $y(\tau)$, but I don't know how to handle $y^n$ since I don't know why the $y^{n+1}=2 \cos (\lambda \tau) y^{n}-y^{n-1}-\frac{\sin (\lambda \tau)}{\lambda} g\left(y^{n}\right)$ is formed as this. Are there any suggestions on how this time integrator formed like above? Or are there any hints on bounding $\left|y\left(t_{n}\right)-y^{n}\right|$? Thanks in advance!
Observations
The integral equation is the variation-of-constants formula for $$ y''+λ^2y+g(y)=0,~~~ y(0)=α,~~ y'(0)=β. $$ A direct finite-difference scheme for that would be $$ \frac{y^{n+1}-2y^n+y^{n-1}}{τ^2}+λ^2y^n=-g(y^n) $$ There is a slight difference in the powers of $τ$ relative to the recursion formula
Applying the integral formula to $t=\pm \tau$ and adding the resulting equations gives \begin{align} y(τ)+y(-τ)&=2\cos(λτ)y(0)-\int_0^τ\frac{\sin(λ(τ-s))}λ g(y(s))\,ds -\int_0^{-τ}\frac{\sin(λ(-τ-s))}λ g(y(s))\,ds \\ &=2\cos(λτ)y(0)-\int_{-τ}^τ\frac{\sin(λ(τ-|s|))}λ g(y(s))\,ds \\ &=2\cos(λτ)y(0)-2g(y(0))\int_0^τ\frac{\sin(λ(τ-s))}λ\,ds \\&~~~~-\int_0^τ\frac{\sin(λ(τ-s))}λ [g(y(s))-2g(y(0))+g(y(-s))]\,ds \\ &=2\cos(λτ)y(0)-2g(y(0))\frac{1-\cos(λτ)}{λ^2}-\int_0^τ\frac{\sin(λ(τ-s))}λ[(g\circ y)''(0)s^2+O(s^4)]\,ds \\ &=2\cos(λτ)y(0)-\frac{4\sin^2(λτ/2)}{λ^2}g(y(0))-\frac{τ^4}{12}(g\circ y)''(0)+O(τ^6) \end{align} One could argue that $4\sin^2(λτ/2)=\sin^2(λτ)+O(τ^4)$, but the square remains, as the coefficient has to be $O(τ^2)$.
Error recursion
The error quantity $e^n=y^n-y(t_n)$ has, with the corrected coefficients, the recursion \begin{align} e^{n+1}-2\cos(λτ)e^n+e^{n-1} &=-\frac{4\sin^2(λτ/2)}{λ^2}[g(y^n)-g(y(t_n))] +\frac{τ^4}{12}(g\circ y)''(t_n)+O(τ^6) \\ &=-\frac{4\sin^2(λτ/2)}{λ^2}g'(y(t_n))e^n+\frac{τ^4}{12}(g\circ y)''(t_n)+O(τ^2|e^n|^2,τ^6) \end{align}
Define additionally $d^n=e^n-e^{n-1}$ so that \begin{align} |e^{n+1}|&\le |e^n|+|d^{n+1}| \\ |d^{n+1}|&\le |d^n|+ \frac{4\sin^2(λτ/2)}{λ^2}\Bigl(λ^2+|g(y(t_n))|\Bigr)|e^n|+\frac{τ^4}{12}|(g\circ y)''(t_n)|+O(τ^2|e^n|^2,τ^6) \\ &\le |d^n|+τ^2L|e^n|+τ^4M+\epsilon^n \end{align} with $L$ and $M$ upper bounds for the factors not depending on the numerical solution, $ϵ^n$ the remainder term. Combining the inequalities gives $$ τ\sqrt{L}|e^{n+1}|+|d^{n+1}|\le (1+τ\sqrt{L}+τ^2L)(τ\sqrt{L}|e^n|+|d^n|) + (1+τ\sqrt{L})(τ^4M+\epsilon^n) $$
Then show inductively that $e^n\in O(τ^2)$, $d^n\in O(τ^3)$, so that $ϵ^n\in O(τ^6)$ and $$ τ\sqrt{L}|e^{n+1}|+|d^{n+1}|\le \begin{aligned}[t]&τ\sqrt{L}(1+τ\sqrt{L}+τ^2L)^n|e^1|\\ &+\frac{τ^3}{\sqrt{L}}\Bigl[(1+τ\sqrt{L}+τ^2L)^n-1\Bigr](M+O(τ^2)) \end{aligned} $$