Question 7.11: EULER’S METHOD FOR SYSTEMS Consider the third-order IVP in E...
EULER’S METHOD FOR SYSTEMS
Consider the third-order IVP in Example 7.9:
3y^{′′′}-y^{′′}-5y^{′}-3y=e^{-x/2} , y(0)=0, y^{′}(0)=-1, y^{′′}(0)=1, 0≤x≤1Using Euler’s method for systems, with step size h = 0.1, find an estimate for y_2 = y(0.2). Confirm by executing the user-defined function EulerODESystem.
In Example 7.9, the IVP was transformed into the standard form of a system of first-order IVPs as
\textbf{u}^{′}=\textbf{f}(x,\textbf{u}), \textbf{u}=\left\{\begin{matrix} u_1 \\ u_2 \\ u_3 \end{matrix} \right\}=\left\{\begin{matrix} y \\ y^{′} \\ y^{′′} \end{matrix} \right\} , \textbf{f}=\left\{\begin{matrix} u_{2} \\ u_{3} \\ \frac{1}{3}\left[u_3+5u_2+3u_1+e^{-x/2}\right] \end{matrix} \right\}, \textbf{u}_0=\left\{\begin{matrix} 0 \\ -1 \\ 0 \end{matrix} \right\}=\left\{\begin{matrix} u_1(0) \\ u_2(0) \\ u_3(0) \end{matrix} \right\}To find y(0.2), we need to find the solution vector \textbf{u}_2 = \textbf{u}(0.2) and then extract its first component, which is y(0.2). By Equation 7.47,
\textbf{u}_{i+1}=\textbf{u}_i+h\textbf{f}(x_i,\textbf{u}_i), i=0,1,2, \cdot \cdot \cdot ,n-1 (7.47)
\textbf{u}_1=\textbf{u}_0+h\textbf{f}(x_0,\textbf{u}_0)
But
\textbf{f}(x_0,\textbf{u}_0)=\left\{\begin{matrix} u_2(x_0) \\ u_3(x_0) \\ \frac{1}{3}\left[u_3(x_0)+5u_2(x_0)+3u_1(x_0)+e^{-x_0/2} \right] \end{matrix} \right\} \overset{x_0=0}{=} \left\{\begin{matrix} -1 \\ 1 \\ \frac{1}{3}(1-5(1)+3(0)+1) \end{matrix} \right\} =\left\{\begin{matrix} -1 \\ 1 \\ -1 \end{matrix} \right\}Therefore,
\textbf{u}_1=\textbf{u}_0+h\textbf{f}(x_0,\textbf{u}_0)=\left\{\begin{matrix} 0 \\ -1 \\ 1 \end{matrix} \right\} +0.1\left\{\begin{matrix} -1 \\ 1 \\ -1 \end{matrix} \right\} =\left\{\begin{matrix} -0.1 \\ -0.9 \\ 0.9 \end{matrix} \right\}In the next step, \textbf{u}_2 = \textbf{u}_1 + h\textbf{f}(x_1, \textbf{u}_1) where
\textbf{f}(x_1,\textbf{u}_1)=\left\{\begin{matrix} u_2(0.1) \\ u_3(0.1) \\ \frac{1}{3}\left[u_3(0.1)+5u_2(0.1)+3u_1(0.1 )+e^{-0.1/2}\right] \end{matrix} \right\} =\left\{\begin{matrix} -0.9 \\ 0.9 \\ \frac{1}{3}(0.9+5(-0.9)+3(-0.1)+e^{-0.1/2} ) \end{matrix} \right\} =\left\{\begin{matrix} -0.9 \\ 0.9 \\ -0.9829 \end{matrix} \right\}Therefore,
\textbf{u}_2=\textbf{u}_1+h\textbf{f}(x_1,\textbf{u}_1)=\left\{\begin{matrix} -0.1 \\ -0.9 \\ 0.9 \end{matrix} \right\}+0.1\left\{\begin{matrix} -0.9 \\ 0.9 \\ -0.9829 \end{matrix} \right\}=\left\{\begin{matrix} \boxed{-0.1900} \\ -0.8100 \\ 0.8017 \end{matrix} \right\}The first component represents y(0.2), thus y(0.2) = −0.19. The results may be confirmed in MATLAB as follows:
>> f = @(x,u)([u(2);u(3);(u(3)+5*u(2)+3*u(1)+exp(-x/2))/3]);
>> x = 0:0.1:1; % 11 mesh points
>> u0 = [0;−1;1];
>> u = EulerODESystem(f,x,u0); % Returns a 3-by-11 matrix
This is a 3 × 11 matrix because there are three state variables and 11 mesh points created. The first row contains the solution estimates for y at the 11 mesh points, the second row y′ and the third row y″. Since we are interested in y(0.2), it is the first row of u that we must retain. In particular,
>> u(1,:)
ans =
Columns 1 through 9
0 −0.1000 [latex]\boxed{−0.1900}[/latex] −0.2710 −0.3440 −0.4099 −0.4698 −0.5245 −0.5751
Columns 10 through 11
−0.6226 −0.6680
The boxed value of y(0.2) agrees with our hand calculations. Knowing the exact value is −0.181324 (truncated), the % relative error associated with our estimate is 4.785%. In order to verify the hand-calculated vectors u_1 and u_2 , we proceed as follows. Since the first column of u represents the initial state vector, we will extract the second and thirds columns to verify the hand calculations:
>> u(:,[2 3])
ans =
−0.1000 −0.1900
−0.9000 −0.8100
0.9000 0.8017
As expected, these exactly match our earlier findings.