7 Nonlinear Equations of Motion

7.1 Limitations of Linear Equations of Motion

Although linearized equations of motion are useful in analyzing controls fixed stability and turning ability of straight line stable ships, they have several limitations. Linearized equations of motion:

cannot be used to analyze the turning ability of directionally unstable ships
cannot accurately predict the behavior of vessels in tight maneuvers
cannot capture speed loss in a turn that can be up to \(40\%\) for high speed destroyers
cannot include coupling effects between speed and steering equations

Nonlinear models are needed to overcome the above limitations. Typically this is done by retaining the higher order terms of the Taylor expansion of hydrodynamic forces that were previously neglected. Popular approaches include terms up to second and third order of the Taylor series expansion.

For most commercial ships, Taylor series beyond third order are generally not considered as they do not contribute significantly to the dynamics in practice. For defense vehicles, it is common to include upto fourth or fifth order terms to capture the intricate dynamics.

7.2 Nonlinear Equations of Motion

The nonlinear speed and steering equations are given by \(\eqref{eq-nonlinear-eom-1}\), \(\eqref{eq-nonlinear-eom-2}\) and \(\eqref{eq-nonlinear-eom-3}\)

\[\begin{align} (m - X_{\dot{u}}) \dot{u} = f_1(u, v, r, \delta) \label{eq-nonlinear-eom-1} \end{align}\]

\[\begin{align} (m - Y_{\dot{v}}) \dot{v} - (Y_{\dot{r}} - m x_G) \dot{r} = f_2(u, v, r, \delta) \label{eq-nonlinear-eom-2} \end{align}\]

\[\begin{align} (I_z - N_{\dot{r}}) \dot{r} - (N_{\dot{v}} - m x_G) \dot{v} = f_3(u, v, r, \delta) \label{eq-nonlinear-eom-3} \end{align}\]

where \(f_1(u, v, r, \delta)\), \(f_2(u, v, r, \delta)\) and \(f_3(u, v, r, \delta)\) are shown in \(\eqref{eq-f1}\), \(\eqref{eq-f2}\) and \(\eqref{eq-f3}\) and contain all the terms upto the third order of the Taylor series expansion.

\[\begin{align} f_1(u, v, r, \delta) =& T - R + X_u (u - U) + X_v v + X_r r + X_{\delta} \delta + X_{uu} (u-U)^2 + X_{uv} (u-U)v \nonumber \\ &+ X_{ur} (u-U)r + X_{u\delta} u \delta + X_{vv} v^2 + X_{vr} vr + X_{v\delta} v \delta + X_{rr} r^2 \nonumber\\ &+ X_{r\delta} r \delta + X_{\delta \delta} \delta^2 + X_{uuu} (u - U)^3 + X_{uuv} (u - U)^2v + X_{uur} (u - U)^2r \nonumber \\ &+ X_{uu\delta} (u - U)^2\delta + X_{uvv} (u - U)v^2 + + X_{uvr} (u - U)vr \nonumber \\ &+ X_{uv\delta} (u - U)v\delta + X_{urr} (u - U)r^2 + X_{ur\delta} (u - U)r\delta \nonumber \\ &+ X_{u\delta \delta} (u - U)\delta^2 + X_{vvv} v^3 + X_{vvr} v^2r + X_{vv\delta} v^2\delta + X_{vrr} vr^2 \nonumber \\ &+ X_{vr\delta} vr\delta + X_{v\delta \delta} v\delta^2 + X_{rrr} r^3 + X_{rr\delta} r^2 \delta + X_{r\delta \delta} r \delta^2 \nonumber \\ &+ X_{\delta \delta \delta} \delta^3 + mvr + mx_G r^2 \label{eq-f1} \end{align}\]

\[\begin{align} f_2(u, v, r, \delta) =& Y_u (u - U) + Y_v v + Y_r r + Y_{\delta} \delta + Y_{uu} (u-U)^2 + Y_{uv} (u-U)v \nonumber \\ &+ Y_{ur} (u-U)r + Y_{u\delta} u \delta + Y_{vv} v^2 + Y_{vr} vr + Y_{v\delta} v \delta + Y_{rr} r^2 \nonumber\\ &+ Y_{r\delta} r \delta + Y_{\delta \delta} \delta^2 + Y_{uuu} (u - U)^3 + Y_{uuv} (u - U)^2v + Y_{uur} (u - U)^2r \nonumber \\ &+ Y_{uu\delta} (u - U)^2\delta + Y_{uvv} (u - U)v^2 + + Y_{uvr} (u - U)vr \nonumber \\ &+ Y_{uv\delta} (u - U)v\delta + Y_{urr} (u - U)r^2 + Y_{ur\delta} (u - U)r\delta \nonumber \\ &+ Y_{u\delta \delta} (u - U)\delta^2 + Y_{vvv} v^3 + Y_{vvr} v^2r + Y_{vv\delta} v^2\delta + Y_{vrr} vr^2 \nonumber \\ &+ Y_{vr\delta} vr\delta + Y_{v\delta \delta} v\delta^2 + Y_{rrr} r^3 + Y_{rr\delta} r^2 \delta + Y_{r\delta \delta} r \delta^2 \nonumber \\ &+ Y_{\delta \delta \delta} \delta^3 - mur + Y_{o} + Y_{ou} (u-U) + Y_{ouu} (u - U)^2 \label{eq-f2} \end{align}\]

\[\begin{align} f_3(u, v, r, \delta) =& N_u (u - U) + N_v v + N_r r + N_{\delta} \delta + N_{uu} (u-U)^2 + N_{uv} (u-U)v \nonumber \\ &+ N_{ur} (u-U)r + N_{u\delta} u \delta + N_{vv} v^2 + N_{vr} vr + N_{v\delta} v \delta + N_{rr} r^2 \nonumber\\ &+ N_{r\delta} r \delta + N_{\delta \delta} \delta^2 + N_{uuu} (u - U)^3 + N_{uuv} (u - U)^2v + N_{uur} (u - U)^2r \nonumber \\ &+ N_{uu\delta} (u - U)^2\delta + N_{uvv} (u - U)v^2 + + N_{uvr} (u - U)vr \nonumber \\ &+ N_{uv\delta} (u - U)v\delta + N_{urr} (u - U)r^2 + N_{ur\delta} (u - U)r\delta \nonumber \\ &+ N_{u\delta \delta} (u - U)\delta^2 + N_{vvv} v^3 + N_{vvr} v^2r + N_{vv\delta} v^2\delta + N_{vrr} vr^2 \nonumber \\ &+ N_{vr\delta} vr\delta + N_{v\delta \delta} v\delta^2 + N_{rrr} r^3 + N_{rr\delta} r^2 \delta + N_{r\delta \delta} r \delta^2 \nonumber \\ &+ N_{\delta \delta \delta} \delta^3 - m x_G ur + N_{o} + N_{ou} (u-U) + N_{ouu} (u - U)^2 \label{eq-f3} \end{align}\]

Note that two subscript term refer to a second derivative term divided by \(2!\) and a three subscript term refers to a third order derivative term divided by \(3!\). The Taylor coefficient \(\frac{1}{n!}\) for the \(n^{th}\) order term is included into the definition of the hydrodynamic coefficient. Thus, \(Y_{uv}\) represents \(\frac{1}{2!}\frac{\partial^2Y}{\partial u \partial v}\) and \(Y_{uvr}\) represents \(\frac{1}{3!}\frac{\partial^3Y}{\partial u \partial v \partial r}\). Note that in the above formulation cross coupling of velocity and acceleration terms are not included. In \(\eqref{eq-f1}\), \(T\) represents the thrust from the propeller and \(R\) denotes the resistance of the hull. \(Y_{o}\), \(Y_{ou}\) and \(Y_{ouu}\) in \(\eqref{eq-f2}\) and \(N_{o}\), \(N_{ou}\) and \(N_{ouu}\) in \(\eqref{eq-f3}\) denote the influence of propeller on sway force and yaw moment respectively.

Exploiting the geometric symmetry of a ship between the port and starboard side, some of the coefficients can be neglected. For example, in the sway equation, it can be observed that \(Y_u = Y_{uu} = Y_{uuu} = 0\) as a change in surge speed alone should not cause any sway force. Similarly, all the coefficients of terms with even powers \(v\), \(r\) and \(\delta\) in sway and yaw equation must be zero. For example, \(Y_{\delta \delta} = 0\) since the sway force due to the rudder must be anti-symmetric about \(\delta = 0\). However, for the surge force all even powers of \(v\), \(r\) and \(\delta\) are retained and all odd powers are neglected. For example, \(X_{vvv} = 0\) as the surge force must develop in the same direction equally for port and starboard sway motions. Neglecting the corresponding hydrodynamic coefficients as discussed here yields the simplified expressions of \(f_1(u, v, r, \delta)\), \(f_2(u, v, r, \delta)\) and \(f_3(u, v, r, \delta)\) as shown in \(\eqref{eq-f1-new}\), \(\eqref{eq-f2-new}\) and \(\eqref{eq-f3-new}\)

\[\begin{align} f_1(u, v, r, \delta) =& T - R + X_u (u - U) + X_{uu} (u-U)^2 \nonumber \\ &+ X_{vv} v^2 + X_{vr} vr + X_{v\delta} v \delta + X_{rr} r^2 \nonumber\\ &+ X_{r\delta} r \delta + X_{\delta \delta} \delta^2 + X_{uuu} (u - U)^3 \nonumber \\ &+ X_{uvv} (u - U)v^2 + X_{uvr} (u - U)vr \nonumber \\ &+ X_{uv\delta} (u - U)v\delta + X_{urr} (u - U)r^2 + X_{ur\delta} (u - U)r\delta \nonumber \\ &+ X_{u\delta \delta} (u - U)\delta^2 + mvr + mx_G r^2 \label{eq-f1-new} \end{align}\]

\[\begin{align} f_2(u, v, r, \delta) =& Y_v v + Y_r r + Y_{\delta} \delta + Y_{uv} (u-U)v \nonumber \\ &+ Y_{ur} (u-U)r + Y_{u\delta} u \delta \nonumber\\ &+ Y_{uuv} (u - U)^2v + Y_{uur} (u - U)^2r \nonumber \\ &+ Y_{uu\delta} (u - U)^2\delta + Y_{vvv} v^3 + Y_{vvr} v^2r + Y_{vv\delta} v^2\delta + Y_{vrr} vr^2 \nonumber \\ &+ Y_{vr\delta} vr\delta + Y_{v\delta \delta} v\delta^2 + Y_{rrr} r^3 + Y_{rr\delta} r^2 \delta + Y_{r\delta \delta} r \delta^2 \nonumber \\ &+ Y_{\delta \delta \delta} \delta^3 - mur + Y_{o} + Y_{ou} (u-U) + Y_{ouu} (u - U)^2 \label{eq-f2-new} \end{align}\]

\[\begin{align} f_3(u, v, r, \delta) =& N_v v + N_r r + N_{\delta} \delta + N_{uv} (u-U)v \nonumber \\ &+ N_{ur} (u-U)r + N_{u\delta} u \delta \nonumber\\ &+ N_{uuv} (u - U)^2v + N_{uur} (u - U)^2r \nonumber \\ &+ N_{uu\delta} (u - U)^2\delta + N_{vvv} v^3 + N_{vvr} v^2r + N_{vv\delta} v^2\delta + N_{vrr} vr^2 \nonumber \\ &+ N_{vr\delta} vr\delta + N_{v\delta \delta} v\delta^2 + N_{rrr} r^3 + N_{rr\delta} r^2 \delta + N_{r\delta \delta} r \delta^2 \nonumber \\ &+ N_{\delta \delta \delta} \delta^3 - m x_G ur + N_{o} + N_{ou} (u-U) + N_{ouu} (u - U)^2 \label{eq-f3-new} \end{align}\]

Notice that the coefficients of the higher order terms in the above equations are known as hydrodynamic coefficients instead of hydrodynamic derivatives (as was done for their linear counterparts). This is because the higher order coefficients are obtained by curve fitting of the data rather than by evaluating slopes. For example, \(Y_{vvv}\) is evaluated by curve fitting the \(Y\) vs \(v\) data rather than evaluating the third order slope of the \(Y\) vs \(v\) curve at \(v=0\).

The nonlinear speed and steering equations along with the kinematic equations can be rearranged to express in the state space form (discussed in detail in later chapters) as shown in \(\eqref{eq-state-space-neom}\).

\[\begin{align} \dot{u} &= \frac{f_1(u, v, r, \delta)}{(m - X_{\dot{u}})} \nonumber \\ \dot{v} &= \frac{(I_z - N_{\dot{r}})f_2(u, v, r, \delta) + (Y_{\dot{r}} - m x_G)f_3(u, v, r, \delta)}{(m - Y_{\dot{v}})(I_z - N_{\dot{r}}) - (N_{\dot{v}} - m x_G)(Y_{\dot{v}} - m x_G)} \nonumber \\ \dot{r} &= \frac{(N_{\dot{v}} - m x_G)f_2(u, v, r, \delta) + (m - Y_{\dot{v}})f_3(u, v, r, \delta)}{(m - Y_{\dot{v}})(I_z - N_{\dot{r}}) - (N_{\dot{v}} - m x_G)(Y_{\dot{v}} - m x_G)} \label{eq-state-space-neom}\\ \dot{x}_{0O} &= u \cos(\psi) - v \sin(\psi) \nonumber \\ \dot{y}_{0O} &= u \sin(\psi) + v \cos(\psi) \nonumber \\ \dot{\psi} &= r \end{align}\]

7.3 Solving Nonlinear Equations of Motion

Together these coupled six equations can be solved to obtain the complete solution of the maneuvering motions in calm water for specified control input \(\delta\). As analytical solutions may not exist in most cases, these equations are usually solved numerically.

In general these equations can be solved using implicit or explicit schemes. Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one. Explicit schemes are fast as the next states are evaluated purely from the knowledge of the current states. Implicit schemes are more computationally expensive as a matrix inversion is usually required to solve for the next states. For schemes achieving the same order of accuracy, implicit schemes are numerically more stable than explicit schemes. Thus, the time step \(h\) can be larger in implicit schemes as compared to explicit schemes.

7.3.1 Euler Scheme

Considering the state space \(\{x\} = [u \quad v \quad r \quad x_{0O} \quad y_{0O} \quad \psi]^T\), the equations of motion depicted in \(\eqref{eq-state-space-neom}\) can be expressed as

\[\begin{align} \{\dot{x}\} = f(t, \{x\}) \end{align}\]

Here, \(\{x\}\) is the unknown function of time. Euler scheme is the simplest explicit scheme where the next state is obtained from the current state by extrapolation along the slope calculated at the current state. Discretizing time \(t_{n+1} = t_{n} + h\) for \(n = 0, 1, 2, ...\), the solution \(\{x\}_{n+1}\) from Euler scheme is given by \(\eqref{eq-euler-scheme}\).

\[\begin{align} \{x\}_{n+1} = \{x\}_{n} + h f\left(t_n, \{x\}_n\right) \label{eq-euler-scheme} \end{align}\]

Given the initial condition \(\{x(t_0)\} = \{x_0\}\), the state space at the next time steps can be evaluated using \(\eqref{eq-euler-scheme}\).

7.3.2 Runge-Kutta \(4^{th}\) order Scheme

One of the popular explicit schemes is the Runge-Kutta \(4^{th}\) order scheme. Discretizing time \(t_{n+1} = t_{n} + h\) for \(n = 0, 1, 2, ...\), the solution \(\{x\}_{n+1}\) is given by

\[\begin{align} \{x\}_{n+1} = \{x\}_{n} + \frac{h}{6} \left(\{k_1\} + 2\{k_2\} + 2\{k_3\} + \{k_4\}\right) \label{eq-rk4} \end{align}\]

where

\[\begin{align} \{k_1\} &= f\left(t_n, \{x\}_n\right) \\ \{k_2\} &= f\left(t_n + \frac{h}{2}, \{x\}_n + \frac{h}{2} \{k_1\}\right) \\ \{k_3\} &= f\left(t_n + \frac{h}{2}, \{x\}_n + \frac{h}{2} \{k_2\}\right) \\ \{k_4\} &= f\left(t_n + h, \{x\}_n + h \{k_3\}\right) \end{align}\]

Thus, given the initial condition \(\{x(t_0)\} = \{x_0\}\), the state space at the next time steps can be evaluated using \(\eqref{eq-rk4}\). This scheme too can be thought of as the extrapolation of current state along a slope to reach the next state. However, unlike the Euler scheme that estimates the slope only at the current state, Runge-Kutta scheme estimates the slope as a weighted average of slopes evaluated at the beginning, midpoint and end of the interval.

There are variants of this scheme that adopt an adaptive time step too to increase the speed (read more here). This adaptive time version is available in most standard programming languages such as MATLAB or Python and can be readily used. In MATLAB, the Runge-Kutta scheme is provided by the function ode45. In Python, the corresponding function is solve_ivp provided through the scipy library.

In the explicit methods as discussed above, the required level of accuracy can be achieved using a sufficiently small time step \(h\). The Euler scheme with a well chosen time step is good enough to simulate the typical maneuvers of a ship. This is due to the fact that the derivatives \(\dot{u}\), \(\dot{v}\) and \(\dot{r}\) typically vary slowly with time due to the large inertia of the ships.

7.4 Memory Effects

Notice that in the Taylor expansion considered for the hydrodynamic forces, a dependence on time was not assumed. This means that the hydrodynamic forces and moments at an instant depend only on the current values of the motion, velocity and acceleration of the vessel and does not retain a memory of the past. Such a process is referred to as a Markov process. Note that in general hydrodynamic forces do not conform to the Markov assumption as waves generated (at the air water interface) at a current time do have effect on its dynamics at a later time. If the motions of the vehicle are sufficiently slow, the memory effects can be disregarded. Since maneuvering motion time scales are much longer than seakeeping motion time scales (occurring due to the presence of waves), the Markov assumption can be made for the maneuvering motions.

Notice that this is also the reason that the frequency dependence is disregarded for the added mass terms (\(X_{\dot{u}}\), \(Y_{\dot{v}}\), \(Y_{\dot{r}}\), \(N_{\dot{v}}\) and \(N_{\dot{r}}\)) and the damping terms (\(X_{u}\), \(Y_{v}\), \(Y_{r}\), \(N_{v}\) and \(N_{r}\)). The frequency dependent added mass terms can be calculated at the zero frequency using frequency domain panel method codes like WAMIT, ANSYS AQWA, Nemoh (open source), HydRA (free to use) etc. Due to this slow motion assumption, the hydrodynamic derivatives are also sometimes referred to as slow derivatives.

7.5 Exercises

Explain why \(Y_u = Y_{uu} = Y_{uuu} = 0\) for a ship with port-starboard symmetry.
Why are higher-order terms in nonlinear maneuvering equations called “hydrodynamic coefficients” rather than “hydrodynamic derivatives”?
For a ship with \(L = 120\) m, if \(Y'_{vvv} = -0.3\), calculate the dimensional coefficient \(Y_{vvv}\). Assume \(\rho = 1025\) kg/m\(^2\).
For a ship with \(X'_{vv} = -0.06\) and \(X'_{rr} = -0.03\), compare the relative magnitudes of non-dimensional surge force due to these hydrodynamic coefficients for pure sway motion \(v' = 0.1\) versus pure yaw motion \(r' = 0.1\).
Explain why even powers of \(v\), \(r\), and \(\delta\) are retained in the surge equation but neglected in sway and yaw equations.
A ship has \(Y'_{vvv} = -0.4\) and \(Y'_{vrr} = -0.2\). For a non-dimensional sway velocity \(v' = 0.15\) and yaw rate \(r' = 0.1\), compare the contributions of these third-order terms to the total sway force.
A ship has the following non-dimensional coefficients:
- \(Y'_v = -0.35\), \(Y'_r = 0.12\), \(Y'_\delta = 0.15\)
- \(Y'_{vvv} = -0.28\), \(Y'_{vrr} = -0.12\)
- \(Y'_{vr\delta} = -0.08\) Calculate the total dimensional sway force in MN when \(v' = 0.2\), \(r' = 0.15\), and \(\delta = 0.1\) rad. You may assume \(L = 200\) m, \(U = 12\) m/s, and \(\rho = 1025\) kg/m\(^3\).
For a ship with \(L = 200\) m and \(U = 12\) m/s, analyze how the ratio of nonlinear to linear contributions in yaw moment changes as rudder angle increases from \(5^{\circ}\) to \(35^{\circ}\) in \(15^{\circ}\) increments. Use:
- \(N'_\delta = 0.15\)
- \(N'_{\delta\delta\delta} = 0.08\)
A ship has both \(Y'_{vvv} = -0.28\) and \(Y'_{vvr} = -0.12\) terms in its sway equation. For a ship with \(L = 200\) m, \(U = 12\) m/s, and \(\rho = 1025\) kg/m\(^3\), calculate the yaw rate \(r\) in rad/s if the sway velocity \(v = 0.2\) m/s and the sway force due to both hydrodynamic coefficients combined is \(Y = 100\) kN.
Using the Euler scheme with time step \(h = 1\) s, calculate the first two time steps for a ship’s yaw rate \(r(t)\) given:
- Initial yaw rate \(r(0) = 0\) rad/s
- \(\dot{r} + 0.02r = 0.05\delta\) where \(\delta = 15^\circ\) (constant)
Compare with the analytical solution at \(t = 2\) s.
For a ship using the RK4 scheme, if the state vector is \(\{x\} = [u \quad v \quad r]^T\), estimate the computational effort (number of function evaluations) needed for:
1. One time step using Euler scheme
2. One time step using RK4 scheme
3. Simulating a 100-second maneuver with \(h = 0.25\) s for both schemes
Discuss the trade-off between accuracy and computational cost.
A ship starting from rest (\(u = 0\) m/s) is given a constant thrust \(T = 10\) kN. Using the surge equation with \(X_u = -50\) kN/(m/s), \(m = 5000\) tonnes, and \(X_{\dot{u}} = -250\) tonnes:
1. Write the differential equation for \(u(t)\)
2. Calculate \(u(0.1)\) using one step of Euler scheme
3. Calculate \(u(0.1)\) using one step of RK4 scheme
4. Compare with the analytical solution
Assume all other forces are negligible and that the nominal operating condition is \(u = 0\) m/s.
Consider a simplified yaw equation \(\dot{r} + 0.5r = 0.1\delta\) with initial condition \(r(0) = 0\) and constant rudder angle \(\delta = 10^\circ\):
1. Write the explicit Euler scheme for this equation
\[\begin{align} r(t_{n+1}) = r(t_n) + h \dot{r}(t_n) \nonumber \end{align}\]
1. Write the implicit Euler scheme for this equation
\[\begin{align} r(t_{n+1}) = r(t_n) + h \dot{r}(t_{n+1}) \nonumber \end{align}\]
1. Calculate \(r(0.5)\) using both schemes with time step \(h = 0.25\) s
2. Compare with the analytical solution and discuss which scheme is more accurate
A ship executes a turning maneuver where its yaw rate follows Nomoto’s model with \(T = 10\) s and \(K = 0.1\) /s. Assuming a constant rudder angle \(\delta = 10^\circ\):
1. Using initial heading \(\psi(0) = r(0) = 0\), calculate \(\psi(2)\) using:
  - Analytical integration
  - Euler scheme with \(h = 1\) s
  - RK4 scheme with \(h = 1\) s
2. Calculate the percentage error in heading angle at \(t = 1\) s for both numerical schemes

Answer Key

A ship with port-starboard symmetry has \(Y_u = Y_{uu} = Y_{uuu} = 0\) because a change in surge speed alone should not cause any sway force.
Higher-order terms in nonlinear maneuvering equations are called ‘hydrodynamic coefficients’ rather than ‘hydrodynamic derivatives’ because they are obtained by curve fitting of the data rather than by evaluating slopes.
\(Y_{vvv} = -2142146.92\) N/(m\(^3\)/s\(^3\))
During pure sway motion, \(X' = -0.0006\) and during pure yaw motion, \(X' = -0.0003\)
Even powers of \(v\), \(r\), and \(\delta\) are retained in the surge equation because the surge force must develop in the same direction equally for port and starboard sway motions. In contrast, they are neglected in sway and yaw equations because the sway force and yaw moment must be anti-symmetric about the nominal operating condition. A positive and negative sway velocity will result in equal and opposite sway forces. Thus the odd powers are retained in the sway and yaw equations.
The contribution of the third-order terms to the total sway force is \(Y_1 = -0.00\) due to \(Y'_{vvv}\) and \(Y_2 = -0.00\) due to \(Y'_{vrr}\)
The total dimensional sway force is \(Y = -118.14\) MN
The ratio of nonlinear to linear contributions in yaw moment changes from \(0.41\%\) at \(5^{\circ}\) to \(6.50\%\) at \(20^{\circ}\) and to \(19.90\%\) at \(35^{\circ}\)
The yaw rate is \(r = -0.0633\) rad/s
The yaw rate at \(t = 2\) s is \(r = 0.0259\) rad/s. The analytical solution is \(r = 0.0257\) rad/s
The computational effort for one time step using Euler scheme is \(1\) function evaluation.

The computational effort for one time step using RK4 scheme is \(4\) function evaluations.

The computational effort for a 100-second maneuver with \(h = 0.25\) s for Euler scheme is \(400\) function evaluations.

The computational effort for a 100-second maneuver with \(h = 0.25\) s for RK4 scheme is \(1600\) function evaluations.

The trade-off between accuracy and computational cost is that RK4 scheme is more accurate but requires more function evaluations.
The differential equation for \(u(t)\) is \((m - X_{\dot{u}})\dot{u} = X_u u + T\)

The surge velocity at \(t = 0.1\) s using Euler scheme is \(u = 0.1905\) m/s

The surge velocity at \(t = 0.1\) s using RK4 scheme is \(u = 0.1904\) m/s

The surge velocity at \(t = 0.1\) s using analytical solution is \(u = 0.1904\) m/s
The yaw rate at \(t = 0.5\) s using explicit Euler scheme is \(r = 0.008181\) rad/s

The yaw rate at \(t = 0.5\) s using implicit Euler scheme is \(r = 0.007326\) rad/s

The yaw rate at \(t = 0.5\) s using analytical solution is \(r = 0.007721\) rad/s

The implicit Euler scheme is more accurate than the explicit Euler scheme
The heading angle at \(t = 2\) s using analytical solution is \(\psi = 0.003269\) rad

The heading angle at \(t = 2\) s using Euler scheme is \(\psi = 0.001745\) rad

The heading angle at \(t = 2\) s using RK4 scheme is \(\psi = 0.003269\) rad

The percentage error in heading angle at \(t = 2\) s using Euler scheme is \(46.61\%\).

The percentage error in heading angle at \(t = 2\) s using RK4 scheme is \(0.00\%\).