# Initial Value Problem for ODE: single-step methods
The solution process of an ODE is performed by steps. In a single-step approach, the solution of the following step depends on the current one. In a multiplestep method, the solution of the next step is calculated from several steps. A multiplestep approach is more accurate, you can think of a similar construct as the higher-order derivatives treated previously.
For first order ODEs with the general form:
$$
\frac{dy}{dx}=f(x,y)
$$
**One initial condition is needed to find only one solution!** Without it, there would be an infinite number of solutions possible. The initial condition is:
$$
y(x_0)=y_0
$$
%% Cell type:markdown id:d1f46802 tags:
# Forward Euler
The simplest initial Value Numerical Integration method is the Forward Euler, which is an explicit method. Although simple, it contains the basic characteristics as more advanced and accurate methods. It looks as:
$$
x_{i+1}=x_i+\Delta x
$$
$$
y_{i+1}=y_i+\Delta x \cdot slope \rvert_i
$$
where $\Delta x$ is the step size and the slope is a constant that approximates the rate of change of $y$ with respect to $x$ (a.k.a. the derivative) in the interval $x_i$ to $x_{i+1}$. The solution starts at $i=0$ given by the initial condition, then $i$ is increased to 1 where the values are calculated using the previous equations. This loop continues until the points cover the desired domain. **The computation of the slope is the key difference between single step methods. For forward euler, the slope is computed with the values at the current step!**
```{figure} figs/explicit_euler.png
:name: explicit_euler
Illustrating the forward Euler method
```
%% Cell type:markdown id:9e6a356f tags:
## Forward Euler Example
Consider the following equation:
$$
\frac{dy}{dx}= y'=-y
$$
We want to know the solution in the domain $[a=0,b=30]$. The **initial condition** is $y(x_0)=1$ and the step size is $\Delta x = 0.2$.
The discretization of the differential equation transforms the problem into an algebraic one. Following the formula of explicit Euler for the discretization yields:
$$
y_{i+1}= y_i + \Delta x *y'_i
$$
In the following code, Forward Euler is implemented.
f(x_{i+1}) = f(x_i) + \Delta x f'(x_i) + \mathcal{O}(\Delta x)
$$
```
:::
%% Cell type:markdown id:9d2bfea3 tags:
## Error analysis
There are two types of errors: round-off and truncation. The **round-off** error is due to the computer's limitation to represent a floating number (decimal). To illustrate it consider the following: the difference between 1 and 0.9 is 0.1. If you subtract 0.1 from itself, you should obtain 0, even if this subtraction is repeated ten thousand times. The following code demonstrates this behavior, using floating-point numbers with varying precision.
%% Cell type:code id:74f6982f tags:
``` python
defaccumulated_error(a,b,solution,iterations):
error=solution-np.abs(a-b)
error_accumulated=0
foriinrange(iterations):
error_accumulated=error_accumulated+error
returnerror_accumulated
print('Error using 16 bits of memory = ',accumulated_error(np.float16(1.),np.float16(0.9),np.float16(0.1),10000))
print('Error using 32 bits of memory = ',accumulated_error(np.float32(1.),np.float32(0.9),np.float32(0.1),10000))
print('Error using 64 bits of memory = ',accumulated_error(np.float64(1.),np.float64(0.9),np.float64(0.1),10000))
```
%% Output
Error using 16 bits of memory = -1.220703125
Error using 32 bits of memory = -0.00022351741790771484
Error using 64 bits of memory = 2.7755575615628914e-13
%% Cell type:markdown id:ee4c7ece tags:
As you can see, this simple operation can give discernible errors. Imagine a computation spanning 100 years with a time step of seconds and more complex operations: the round-off error will be present! As you can see, this can be reduced by increasing the precision or the number of digits used to represent numbers but be careful, this is not free as the memory the computer uses increases as well as the computation time.
The truncation error is related to the method chosen for the slopes approximation. This can be obtained using our reliable TSE as this gives the exact solution. Therefore **the truncation error per step** is
The error that is introduced in each step of the numerical solution ideally does not to increase as the solution advances. Under a well posed and proper solution, the error is expected to reduce with smaller steps. In some cases, the error increases without bound as the solution advances (even with smaller steps): **the solution becomes unstable**. The stability depends on the numerical method, the step size and the behavior of the differential equation. Therefore, the stability conditions will differ when applying the same numerical method to different equations.
Let's consider the **stability for Forward Euler** applied to a more general form of the problem above.
$$
\frac{dy}{dx} = y'=-\alpha y
$$
With initial condition $y(0)=1$ and $\alpha>0$, the exact solution is:
Comparing this last equation with the exact solution, it can be seen that the term $(1-\alpha \Delta x)^n$ in the numerical solution is approximating the term $e^{-\alpha \Delta x_i}$ in the exact solution. The latter tends to decay with larger values of $x_i$ and positive $\alpha$. To make sure that the term $(1-\alpha \Delta x)^n$ decays with larger $n$ values $1-\alpha \Delta x$ should be less than $|1|$, i.e.,
```{admonition} Click here for a more detailed derivation
:class: tip, dropdown
$$
-1 <1-\alpha \Delta x < 1
$$
Subtract 1 of both sides
$$
-2 <-\alpha \Delta x < 0
$$
Divide by -1
$$
0 <\alpha \Delta x < 2
$$
```
This is the stability criterion. If not complied, then the solution will be unstable. It is said that **Explicit Euler is conditionally stable! This is true for every explicit numerical method.**
%% Cell type:markdown id:d76d9cd7 tags:
:::{card} Exercise
Lets go back to the last example described in the code above.
$$
\frac{dy}{dx}= y'=-y
$$
We want to know the solution in the domain $[a=0,b=30]$. The **initial condition** is $y(x_0)=1$ and the step size is $\Delta x = 0.2$.
What is the step size at which the problem described using the Forward Euler becomes unstable? Feel free to use the interactive figure below and change dx to test what happens when the step size becomes larger.
Click `rocket` -->`Live Code` to interact with the plot below
You can see in the interactive figure above that when the stepsize becomes larger the forward euler solution becomes less accurate. You can also see that when the step size becomes bigger then 2 the Forward Euler solution starts to blow up. This is because of the stability criteria of the Forward Euler.
$$
0 <\alpha \Delta x < 2
$$
When $\alpha=1$, then $\Delta x <2$(Notethat$\Deltax$cannotbenegativeandthereforewillalwaysbe> 0 )
```
%% Cell type:markdown id:20307b08 tags:
# Backward Euler
%% Cell type:markdown id:aad207a0 tags:
Just as Forward Euler is the simplest Explicit method, **Backward Euler method is the simplest implicit method**. They look very similar but the slope is computed at the next step! That's where the implicit term comes from: the slope depends on the unknown value.
$$
x_{i+1}=x_i+\Delta x
$$
$$
y_{i+1}=y_i+\Delta x \cdot slope \rvert_{i+1}=y_i+\Delta x \cdot y'_{i+1}
$$
Let's consider the same problem as before:
$$
\frac{dy}{dx}= y'=-\alpha y
$$
We want to know the solution in the domain $[a=0,b=30]$. The **initial condition** is $y(x_0)=1$ and the step size is $\Delta x = 0.2$.
The discretization of the differential equation transforms the problem into an algebraic one. Following the formula of implicit Euler, the discretization yields:
$$
y_{i+1}= y_i - \Delta x \cdot \alpha y_{i+1}
$$
then
$$
y_{i+1}+ \Delta x \cdot \alpha y_{i+1} = y_i
$$
finally
$$
y_{i+1}= y_i/(1+\Delta x \cdot \alpha)
$$
In the following code, Implicit Euler is implemented.
You can see that the result is similar to Forward Euler, except that the curve for Backward Euler is slightly above than the exact solution. This makes sense as the derivative is taken to be in front, so the slope is initially underestimated while the explicit solution initially overestimated the slope.
The round-off errors and truncation error are similar. Backward Euler is also a first-order method. But what about the stability?
:::
%% Cell type:markdown id:fda50b29 tags:
:::{card} Exercise:
Modify the step size in the Implicit Euler code, try to make the solution unstable. What do you notice?
This last term to the power $n$ approximates the decaying exponential in the exact solution, just like in the Explicit Euler method. Here, the only condition that must be satisfied for stability is
$$\left|\frac{1}{1+\alpha \Delta x}\right| < 1$$
This can be expressed as:
$$ -1 <\frac{1}{1+\alpha \Delta x} < 1$$
Lets first look at the right side condition $\frac{1}{1+\alpha \Delta x} < 1$.
If we multiply by the term $1+\alpha \Delta x$, we get:
$$
1 < 1+\alpha \Delta x
$$
As $\alpha$ is positive and the step size must be larger than 0, the right side of the criterion is always true.
Lets now look at the left side condition $-1 <\frac{1}{1+\alpha \Delta x}$:
We again multiply by $1+\alpha \Delta x$:
$$
-1 - \alpha \Delta x < -1
$$
Again this left side criterion always holds.
Hence, **an implicit Backward Euler scheme is unconditionally stable.** This is also the case for more advanced implicit schemes.
%% Cell type:markdown id:560fbf20 tags:
:::{card}
:::{card} Exercise
Use the forward Euler to approximate the solution to the initial-value problem.
