.. _ch02-2.7: ======================================================== [7] Numerical Approximations: Euler's Method ======================================================== * In this section, we learn how to solve DEs computationally, or in other words, by writing a computer code to solve DEs numerically. * Many practical applications using mathematical modelings with DEs (ODEs and PDEs) often do not have analytical solutions. This means that we do not know how to solve those problems in practice. * The methods and examples we learned so far are indeed only a fractional part of DEs in which simple analytical ways are available to solve them. * When analytical theories/methods are not avialble, we then need to rely on computer models to study such general, and more challenging DEs. * One easy way, due to `Euler `_ in about 1768, is based on using tangent lines of the DEs, and the method is called the Euler method (or the tangent line method). To be precise, the method we are going to learn is called the *Forward* Euler Method in order to represent the discrete approximation of the method is *forward* in time. That is to say, by considering the definition of derivatives, we see that, for :math:`y'=f(t,y)`, .. math:: \lim_{h\to 0} \frac{y(t+h) - y(t)}{h} = \frac{dy}{dt} = f(t,y). * If we choose :math:`h` very small then we have .. math:: :label: forwardEuler \frac{y(t+h) - y(t)}{h} \approx \frac{dy}{dt} = f(t,y). * This is called the *forward differencing* approximation of :math:`\frac{dy}{dt}`, whereas the *backward differencing* approximation can be similarly written as .. math:: \frac{y(t) - y(t-h)}{h} \approx \frac{dy}{dt} = f(t,y). .. admonition:: The (forward) Euler's Method for ODEs Given the IVP: .. math:: :nowrap: \begin{equation} \left \{ \begin{array}{ccc} && \frac{dy}{dt} = f(t,y)\nonumber \\ && y(t_0) = y_0, \end{array} \right. \end{equation} The Euler method takes the following steps: 1. Choose the initial point :math:`(t_0,y_0)` of the IVP. 2. Introduce a discrete notation :math:`t_1` to write :math:`t_0 + h`, and in general, .. math:: t_n = t_0 + n h. 3. Introduce another discrete notation :math:`y_n` to write .. math:: y_n = y(t_n) = y(t_0 + n h). 4. Then Eq. :eq:`forwardEuler` can be written as a discrete difference equation: .. math:: \frac{y_{1} - y_0}{h} \approx f(t_0,y_0), or, equivalently, .. math:: y_1 \approx y_0 + h f(t_0,y_0) 5. Assume :math:`h` is chosen very small enough so that we approximately write .. math:: y_1 = y_0 + h f(t_0,y_0), and in the similar way, we successively obtain the :math:`(n+1)` th term, .. math:: y_{n+1} = y_{n} + h f(t_{n},y_n). .. figure:: ./_figures/fig2.7.2.png :align: center :scale: 50% Examples 1 and 2: ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Solve the IVP .. math:: \frac{dy}{dt} = 3 - 2t - 0.5y, y(0) = 1, using the Forward Euler method. Note that using the integrating factor :math:`e^{0.5t}` we can obtain the exact solution .. math:: y = \phi(t) = 14 - 4t - 13 e^{-0.5t}. **An example Matlab code:** Click to download the sample code :download:`ForwardEuler.m <./code/ForwardEuler.m>` .. literalinclude:: ./code/ForwardEuler.m :language: matlab :linenos: **Questions** 1. Study and understand the code, and run it. 2. Modify the time step size :math:`h` by increasing or decreasing within the stability range of .. math:: 0 \le h \le 4. What happens when you increase/decrease your timestep :math:`h`? 3. Modify the time step size :math:`h` by increasing more than the stability range so that .. math:: h > 4. What happens when your timestep :math:`h` is bigger than the stability limit 4? 4. The current code implements to exit the temporal evolution based on the maximum timestep, e.g., ``n_step``. Modify the exit condition so that the code stops when ``t=t_max``, where ``t_max`` is any arbitrary time which is not necessarily a multiple of ``n_step`` as in the current implementation. (hint: use ``while`` loop instead of using ``for`` loop.) 5. Check if you can reproduce the numbers in Table 2.7.1 and Table 2.7.2. .. figure:: ./_figures/tab2.7.1.png :align: center :scale: 50% .. figure:: ./_figures/fig2.7.3.png :align: center :scale: 50% .. figure:: ./_figures/tab2.7.2.png :align: center :scale: 50% .. note:: * We do not learn how to obtain a stability limit for timestep :math:`h` in this class. Such topics are covered in more advanced numerical ODE/PDE classes (e.g., `AMS 213B `_). * For this problem, as long as you keep the time step size :math:`h` to satisfy the stability bound, .. math:: 0 \le h \le 4, the numerical solutions of the Forward Euler method are stable. Example 3: ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Solve the IVP .. math:: \frac{dy}{dt} = 4 - t + 2y, y(0) = 1, using the forward Euler method. **Questions** 1. Find the exact solution. 2. Modify the code in Examples 1 and 2 to solve this problem. 3. Repeat the steps in Example 1 and 2 and see if the solutions converge by visually inspecting the departure between the numerical and exact solutions, in particular, for large :math:`t`. .. note:: * For this problem, the mathematical stability bound of the time step size :math:`h` becomes .. math:: -1 \le h \le 0, which implies the Forward Euler method **always fails** to be stable, no matter how :math:`h>0` is small, for the given ODE. * The unstably divergent solution behaviors are shown in Table 2.7.3, in which the differences between the exact and numerical values at times :math:`t=0,1,2,3,4,5` are much bigger than Examples 1 and 2, for all choices of :math:`h`. * In this situation, one has to use a different numerical method, such as the backward Euler method (harder than using the forward Euler method in general), in order to stably solve the problem numerically. .. figure:: ./_figures/tab2.7.3.png :align: center :scale: 50% .. figure:: ./_figures/fig2.7.4.png :align: center :scale: 50%