How do I determine the order of consistency of a method?
I do see a lot of approaches with the Taylor expansion, but I don't really understand them and I don't really see the Taylor Expansion in them.
For example:
Given is an ODE $ y' = f(t,y) $ with $ y(0)=y_0 $ and the implicit Euler method to solve the ODE as $ y_{j+1} = y_j +hf(t_{j+1}, y_{j+1}) $, where $ y_j=y(t_j) $ and $ y_{j+1}=y(t_{j+1}) $.
When I want to determine the order of consistency I would start like this:
$$(1)h\tau_{j+1} = y(t_{j+1})-y_h(t_{j+1}) $$
$$(2)=y(t_{j+1})-y(t_j)-hf(t_{j+1}, y_h(t_{j+1}) $$
$$(3)=y(t_{j+1})-y(t_j)-hf(t_{j+1}, y(t_{j+1}) – h\tau_{j+1})$$
Accoring to a solution this would be the next step made with Taylor Expansion:
$$(4)=y(t_{j+1})-y(t_j)-hf(t_{j+1}, y(t_{j+1})) + h\frac{\partial f}{\partial y}(t_{j+1},\xi)h\tau_{j+1}$$
So my question here is: How do I come exactly from (3) to (4)? How do I use the Taylor Expansion correctly here to get the correct result?
Best Answer
In reality, the goal is estimate the global error, i.e., the difference $$e_j = y(t_j) - v_j$$ between the exact solution $y(t_j)$ and the corresponding approximation $v_j$ which has been computed using, say, Euler's implicit method.
To that end, we seek an inequality which is satisfied by the error. We have the ordinary differential equation $$y'(t) = f(t,y(t))$$ and we have the discrete equation, in this case $$v_{j+1} = v_j + h f(t_{j+1}, v_{j+1}).$$ While $y_j = y(t_j)$ is rarely identical to $v_j$ we hit upon the idea to investigate the extent to which $y_j$ satisfies Euler's scheme. By Taylor's theorem we have $$y_j = y(t_j) = y(t_{j+1} - h) = y(t_{j+1}) - y'(t_{j+1})h + \frac{1}{2}y''(\xi_j) h^2, $$ where $\xi_j$ is between $t_j$ and $t_{j+1}$. It follows that $$ y(t_{j+1}) = y(t_j) + h f(t_{j+1}, y(t_{j+1}))h - \frac{1}{2}y''(\xi_j)h^2$$ or equivalently $$ y_{j+1} = y_j + h f(t_{j+1}, y_{j+1}) - \frac{1}{2}y''(\xi_j)h^2.$$ This is very nearly Euler's implicit method. At this point, the order of consistency can be read off. Some authors prefer to subtract $1$ from $2$, the order the term $- \frac{1}{2}y''(\xi_j)h^2$. Others prefer to work with an equivalent equation, namely $$ \frac{y_{j+1}-y_j}{h} - f(t_{j+1}, y_{j+1}) = - \frac{1}{2}y''(\xi_j)h $$ and simply read off the order of the right hand side term. In either case, the order of consistency of Euler's implicit method is $1$.
We now return to the analysis of the global error. By combining the two equations, we have \begin{align} e_{j+1} &= y_{j+1} - v_{j+1} \\ &= y_j + h f(t_{j+1}, y_{j+1}) - \frac{1}{2}y''(\xi_j)h^2 - (v_j - h f(t_{j+1}, v_{j+1})) \\ &= y_j - v_j + h (f(t_{j+1}, y_{j+1}) - f(t_{j+1}, v_{j+1})) - \frac{1}{2}y''(\xi_j)h^2\end{align} This is the point where it gets exciting. We usually have imposed a Lipschitz condition of the form $$| f(t,y_1) - f(t,y_2) | \leq L |y_1 - y_2 |$$ for all relevant $t$, $y_1$ and $y_2$. This is a very mild requirement which is instrumental in securing the existence and uniqueness of the true solution of the initial value problem under investigation. Here the same condition can be used to bound the error. We have $$|e_{j+1}| \leq (1 + Lh)|e_j| + \frac{1}{2} |y''(\xi_j)| h^2.$$ If we limit ourselves to a closed and bounded interal, say $t \in [0,T]$, then we can bound the absolute value of $y''$ by some constant $M$ and arrive at $$|e_{j+1}| \leq (1 + Lh)|e_j| + \frac{1}{2} M h^2.$$ This can be iterated to yield $$ |e_j| \leq \frac{1}{2} M \left( \frac{e^{LT}-1}{L}\right) h, \quad jh \leq T$$ It is no coincidence that the order of the global error is also $1$.