User:Egm6341.s11.team2.finalreportHW7

=Problem 7.0.1: Linear State Space Model =

Given: The Linear State Space Model (LSSM)
Where LSSM has the general form
 * {| style="width:100%" border="0"

$$ \displaystyle \mathbf{x_{k+1}}_{\color{red}(nx1)} = \mathbf{F}_{\color{red}(nxn)} \mathbf{x_{k}}_{\color{red}(nx1)} $$     (7.0.1)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

And $$ \displaystyle \mathbf{F}_{\color{red}(nxn)} $$ is defined as
 * {| style="width:100%" border="0"

$$ \displaystyle \mathbf{F}_{\color{red}(nxn)} = \mathbf{I}_{\color{red}(nxn)} + \Delta \mathbf{A}_{\color{red}(nxn)} $$     (7.0.2)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

If we choose $$ \displaystyle {\color{red} n=2} $$ and define the matrices in (2.2) as
 * {| style="width:100%" border="0"

$$ \displaystyle \mathbf{I}_{\color{red}(2x2)} = \left[ \begin{matrix} 1 & 0 \\   0 & 1  \\ \end{matrix} \right], \quad \mathbf{\Delta} = 0.02, \quad \mathbf{A}_{\color{red}(2x2)} = \left[ \begin{matrix} -0.2 & 1    \\   -1   & -0.2  \\ \end{matrix} \right] $$     (7.0.3)
 * style="width:95%" |
 * style="width:95%" |
 * 
 * }

In addition we will consider the k-th step of the system $$ \displaystyle \mathbf{x_k} $$ and initial point $$ \displaystyle \mathbf{x_0} $$, in (7.0.1) as
 * {| style="width:100%" border="0"

$$ \displaystyle \mathbf{x_k}_{\color{red}(2x1)} = \left\{ \begin{matrix} x_k^ \\ x_k^ \\ \end{matrix} \right\} \quad and \quad
 * style="width:95%" |
 * style="width:95%" |

\mathbf{x_0}_{\color{red}(2x1)} = \left\{ \begin{matrix} x_k^ \\ x_k^ \\ \end{matrix} \right\} = \left\{ \begin{matrix} 3 \\   -2  \\ \end{matrix} \right\} $$     (7.0.4)
 * 
 * }

1. Run LSSM With Cauchy Random Noise:
====a) Let $$ \displaystyle \mathbf{G} = \alpha \cdot \left\{ \begin{matrix}  1  \\   1  \\ \end{matrix} \right\}_{\color{red}2x1}  $$.====

Hint:
Find a Matlab command to generate $$ \displaystyle \{ \mathbf{\theta_{j}}, j=0,1,2,... \} $$ in single-slit diffraction experiment.

1. Run LSSM With Cauchy Random Noise:
====a) Let $$ \displaystyle \mathbf{G} = \alpha \cdot \left\{ \begin{matrix}  1  \\   1  \\ \end{matrix} \right\}_{\color{red}2x1}  $$.====

c) Using the Matlab command rand to generate $$ \displaystyle \{ \mathbf{w_{j}}, j=0,1,2,... \} $$.
Matlab Code:

d) Plot $$ \displaystyle \{ \mathbf{x_{j}}, j=0,1,2,... \} $$ for $$ \displaystyle \alpha=0.5,1,2 $$.
Shown Below is the Linear State Space Model With Cauchy Noise and $$ \displaystyle {\color{red}\alpha = 0.5} $$:
 * HW701cauchyNoiseAlpha_05.png

Shown Below is the Linear State Space Model With Cauchy Noise and $$ \displaystyle {\color{red}\alpha = 1} $$:
 * HW702cauchyNoiseAlpha_1.png

Shown Below is the Linear State Space Model With Cauchy Noise and $$ \displaystyle {\color{red}\alpha = 2} $$:
 * HW703cauchyNoiseAlpha_2.png

Hint:
Find a Matlab command to generate $$ \displaystyle \{ \mathbf{\theta_{j}}, j=0,1,2,... \} $$ in single-slit diffraction experiment.

=Problem 7.0.2: Mass-Spring-Damper with Gaussian noise and Cauchy noise=

Given: The Following Spring System
An ideal mass-spring-damper system with mass m, spring constant k and viscous damper of damping coefficient c is subject to an oscillatory force u. This can be illustrated in the following image.



Find:
1) Derive equations of motion in terms of d, c, k, m, u.

2) Let $$ {\mathbf{x}} = \left\{ {\begin{array}{ccccccccccccccc} d \\  {\dot d} \end{array}} \right\} = \left\{ {\begin{array}{ccccccccccccccc}   \\ \end{array}} \right\}$$. Find $$\left( {{\mathbf{F}},{\mathbf{G}}} \right) $$

3) Find $${c_{cr}}$$ in terms of k, m st. this system is critically damped.

4) Let k=1, m=1/2, $$ {{\mathbf{x}}_0} = {\left[ {0.8, - 0.4} \right]^T} $$

a) For u=0, plot $${{\mathbf{x}}_k}$$ for $$c = \frac{1}{2}{c_{cr}},{c_{cr}},\frac{3}{2}{c_{cr}}$$.

b) For u=0.5 Gaussian noise and $$c = \frac{3}{2}{c_{cr}}$$, plot $${{\mathbf{x}}_k}$$.

c) For u=0.5 Cauchy noise and $$c = \frac{3}{2}{c_{cr}}$$, plot $${{\mathbf{x}}_k}$$.

Solution:
1) Since the system is subject to an oscillatory force generated from the spring, we have this oscillatory force

and also a damping force from the damper:

Since we can use Newton's Second Law

Therefore we can derive equations of motion as following:

2) From Eq. 5.4.4, we can solve for $${\ddot d}$$ as following:

Therefore we can rewrite $${\dot x}$$ in matrix form

If we make a discretization as

and use forward Euler method

We can have

3) Since we can rewrite Eq. 5.4.4 as

where $$\omega_0 = \sqrt{ k \over m } $$, $$\zeta = { c \over 2 \sqrt{m k} }$$

When $$\zeta = 1$$, the system is said to be critically damped. A critically damped system converges to zero faster than any other and without oscillating. Therefore we have $${c_cr}$$ indicated in Eq. 5.4.12 and can look for damping for more explanations.

4) a) The following matlab codes are used to solve the problem without any noise, and the damping coefficients can be chose to create different cases with under-damping, critical damping or over-damping.



Comment:

Different c values cause different damping cases. From the figure above, we can find out a critically damped system converges to zero faster than any other and without oscillating.

b) The following matlab codes are used to solve the problem with the Gaussian Random Noise and over-damping case.



Comment:

When the system has Gaussian noise, the convergence can be archived, however the final result is not exactly zero since there is still noise all the time and the final result oscillated around zero.

c) The following matlab codes are used to solve the problem with the Cauchy Noise and over-damping case.



Comment:

When the system has Cauchy noise, the same result was generated as Gaussian noise. However the oscillation around the zero is much larger than the case with Gaussian noise.

=Problem 7.1.1: Solve the Logistic Growth ODE Analytically=

Given: The Verhulst Model for Logistic Population Growth
The Verhulst Model, as it pertains to logistic population growth, describes the following system dynamics.
 * LHS: The time derivative $$ \displaystyle \dot{x}(t)= \frac{dx}{dt} $$, is the rate that the population changes with respect to time.
 * RHS: Is directly proportionate to the following:
 * $$ \displaystyle x(t):= $$ the population size at time $$ \displaystyle t$$.
 * $$ \displaystyle \left(1-\frac{x(t)}{x_{max}} \right):= $$ the remaining available resources.
 * $$ \displaystyle x_{max}:= $$ the carrying capacity which is maximum number of inhabitants for an isolated system.
 * $$ \displaystyle r:= $$ the constant of proportionality which serves as the rate constant. We are considering the case of population growth where $$ \displaystyle r>1 $$.

Find: An Analytical Solution for the Population Size at Time $$ t $$

 * Obtain the analytical solution by solving the ODE shown in equation (7.1.1).
 * Use the initial condition $$ \displaystyle x_0 = 3 < \frac{1}{2}x_{max} $$.
 * Use the initial condition $$ \displaystyle x_0 = 9 < \frac{1}{2}x_{max} $$.

Step 1: Separation of Variables
Separating the independent and dependent variables in equation (7.1.1) results in the following expression that can then integrated.

Step 2: Partial Fraction Decomposition
Next we will consider the integrand and decompose its fractional form into two the sum of two elementary integrals.

Multiplying both the LHS and RHS by the denominator on the LHS results in the following, The final expression in (7.1.4) can be written in terms of its respective powers of $$ \displaystyle x $$, forming linear system of two equations with two unknowns. The unknowns we aim to obtain are $$ \displaystyle a $$ and $$ \displaystyle b $$. Since $$ \displaystyle a $$ is coefficient to $$ \displaystyle x^0=1 $$, we immediately obtain,

Next we observe the coefficients of order one $$ \displaystyle x $$ and using our solution from (7.1.5) $$ \displaystyle a=1 $$ to get,

Now that we have a solution to $$ \displaystyle a $$ and $$ \displaystyle b $$ we can substitute them into (7.1.3) to obtain an integrand that is integrable by elementary techniques.

Step 3: Direct Integration
Next we will simply evaluate the integral and perform the algebra necessary in obtaining and explicit expression for $$ \displaystyle x(t) $$ Which gives us the Verhulst Model for logistic population growth.

Step 4: Using the Data Set Parameters for $$ t \in [0,20] $$

 * Population carrying capacity $$ \displaystyle x_{max}=15 $$.
 * Population growth rate $$ \displaystyle r = 1.4 $$.
 * Initial condition $$ \displaystyle x_0^{(1)} = 3 < \frac{1}{2}x_{max} $$.


 * Population carrying capacity $$ \displaystyle x_{max}=15 $$.
 * Population growth rate $$ \displaystyle r = 1.4 $$.
 * Initial condition $$ \displaystyle x_0^{(2)} = 9 > \frac{1}{2}x_{max} $$.

= Problem 7.1.2: Use Hermite-Simpson Algorithm to Numerically Integrate the Verhulst Equation=

Given: The Verhulst Model for Logistic Population Growth
Verhulst(1838) equation

The Verhulst Model, as it pertains to logistic population growth, describes the following system dynamics.
 * LHS: The time derivative $$ \displaystyle \dot{x}(t)= \frac{dx}{dt} $$, is the rate that the population changes with respect to time.
 * RHS: Is directly proportionate to the following:
 * $$ \displaystyle x(t):= $$ the population size at time $$ \displaystyle t$$.
 * $$ \displaystyle \left(1-\frac{x(t)}{x_{max}} \right):= $$ the remaining available resources.
 * $$ \displaystyle x_{max}:= $$ the carrying capacity which is maximum number of inhabitants for an isolated system.
 * $$ \displaystyle r:= $$ the constant of proportionality which serves as the rate constant. We are considering the case of population growth where $$ \displaystyle r>1 $$.

2. Using the initial condition $$ \displaystyle x_0 = 9 > \frac{1}{2}x_{max} $$

 * Obtain the numerical solution by solving the ODE shown in equation (7.2.1.1) for the following parameters.
 * Given the population carrying capacity $$ \displaystyle x_{max}=15 $$.
 * Given the population growth rate $$ \displaystyle r = 1.4 $$.

Solution: Numerical Integration
The following Matlab codes are used to solve Eq. 7.1.2.1

1)$$ x_0 = 3$$



2)$$ x_0 = 9$$



3)comparison between the above cases



Comment
Form the above two cases, we find that these two cases, with different initial condition, can converge to the maximum value exactly. It's interesting that it seems that both cases achieve the maximum at the same time. However, the increase rate is different for they have various initial start point.

=Problem 7.1.3: FE and BE Algorithm Applied to the Integrations of the Logistic Growth Equation=

From the lecture slide Mtg 40

 This problem was solved with referring to S10 homework: S10 Team4's Problem 7.5 

Given
The Verhulst Model as stated in the Problem 7.1.2,

Two cases of initial condition,


 * Case 1: initial condition $$ \displaystyle x_0 = 3 < \frac{1}{2}x_{max} $$.
 * Case 2: initial condition $$ \displaystyle x_0 = 9 < \frac{1}{2}x_{max} $$.

Find
Integrate the logistic equation with increasing step size $$\displaystyle h=\hat{h}\cdot {{2}^{k}}$$ by

(1) Forward Euler algorithm

(2) Backward Euler algorithm

Integrate using Backward Euler Algorithm
Do a scaling such that $$\displaystyle \bar{x}(t)=\frac{x(t)}$$, then the Verhulst Model turns into,

Then by Forward Euler Algorithm,

Case 1: $$\displaystyle {{x}_{0}}=3$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$



Result using Forward Euler Algorithm when $$\displaystyle k=4$$



Result using Forward Euler Algorithm when $$\displaystyle k=7$$



Case 2: $$\displaystyle {{x}_{0}}=9$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$



Result using Forward Euler Algorithm when $$\displaystyle k=4$$



Result using Forward Euler Algorithm when $$\displaystyle k=7$$



Integrate using Backward Euler Algorithm
From Eq 1.3.2 and Backward Euler Algorithm we have,

Case 1: $$\displaystyle {{x}_{0}}=3$$
Result using Forward Euler Algorithm when $$\displaystyle k=1$$



Result using Forward Euler Algorithm when $$\displaystyle k=7$$



Case 2: $$\displaystyle {{x}_{0}}=9$$
Result using Forward Euler Algorithm when $$\displaystyle k=4$$



Result using Forward Euler Algorithm when $$\displaystyle k=7$$



Comment
As we can see in the Foward Euler Algorithm the result become indented when $$\displaystyle k=7$$ reaches 7. This means the Foward Euler Algorithm's solution may become unstable with step size increasing. While the result by Backward Euler Algorithm still remain stable with $$\displaystyle k=7$$.

=Problem 7.2: Using Hermite-Simpson and Newton-Raphson Algorithm to Solve the Bunt Maneuver Problem=

Given
The given equations of motion for the aircraft are


 * {| style="width:100%" border="0" align="left"

\mathbf{\dot{z}} = \begin{Bmatrix} \dot{\gamma} \\ \dot{v} \\ \dot{x} \\ \dot{y} \end{Bmatrix} = \begin{Bmatrix} \frac{T-D}{mv}\sin\alpha + \frac{L}{mv}\cos\alpha - \frac{g}{v}\cos\gamma \\ \frac{T-D}{m}\cos\alpha + \frac{L}{m}\sin\alpha - g\sin\gamma \\ v\cos\gamma \\ v\sin\gamma \end{Bmatrix} = \begin{Bmatrix} f_{1}(\mathbf{z}) \\ f_{2}(\mathbf{z}) \\ f_{3}(\mathbf{z}) \\ f_{4}(\mathbf{z}) \end{Bmatrix} = \mathbf{f}(\mathbf{z}) $$ $$  Here, $$\displaystyle D$$ and $$\displaystyle L$$ are
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\begin{cases} D(y,v,\alpha) = \frac{1}{2}C_{d}\rho v^{2}S_{ref} \\ L(y,v,\alpha) = \frac{1}{2}C_{l}\rho v^{2}S_{ref} \end{cases} $$  where
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }
 * {| style="width:100%" border="0" align="left"

\begin{cases} C_{d}(\alpha) = A_{1}\alpha^2 + A_{2}\alpha + A_{3} \\ C_{l}(\alpha) = B_{1}\alpha + B_{2} \\ \rho(y) = C_{1}y^2 + C_{2}y + C_{3} \end{cases} $$  
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

All required parameters can be obtained from the Subchan and Zbikowski paper. (Subchan and Zbikowski, Optim. Control Appl. Meth. 2007; 28:311-353)
 * {| style="width:100%" border="0" align="left"

\begin{cases} m = 1005 \; kg \\ g = 9.81 \; m/s^2 \\ S_{ref} = 0.3376 \; m^2 \\ A_{1} = -1.9431 \\ A_{2} = -0.1499 \\ A_{3} = 0.2359 \\ B_{1} = 21.9 \\ B_{2} = 0 \\ C_{1} = 3.312 \times 10^{-9} \; kg/m^5 \\ C_{2} = -1.142 \times 10^{-4} \; kg/m^4 \\ C_{3} = 1.224 \; kg/m^3 \end{cases} $$   
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }

Find
Solve the equations of motion of aircraft when the below control inputs $$\displaystyle T(t)$$ and $$\displaystyle \alpha_1(t)$$ and $$\displaystyle \alpha_2(t)$$ are given.


 * {| style="width:100%" border="0" align="left"

T(t) = \begin{cases} 6000 \; N, & \mbox{for }0 \leq t < 27 \; \mbox{and} \; 33 \leq t \leq 40 \\ 1000 \; N, & \mbox{for }27 \leq t < 33 \end{cases} $$ $$ 
 * $$\displaystyle
 * $$\displaystyle
 * $$\displaystyle
 * }
 * }


 * {| style="width:100%" border="0" align="left"

\alpha_1(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.09}{27-21}(t-21) + 0.09 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-21}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ \frac{-0.13+0.2}{40-33}(t-33) - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$ $$  
 * $$\displaystyle
 * $$\displaystyle
 * <p style="text-align:right;">$$\displaystyle
 * }
 * }


 * {| style="width:100%" border="0" align="left"

\alpha_2(t) = \begin{cases} 0.03 \; rad, & \mbox{for }0 \leq t < 21 \\ \frac{0.13-0.09}{27-21}(t-21) + 0.09 \; rad, & \mbox{for }21 \leq t < 27 \\ \frac{-0.2+0.13}{33-21}(t-27) - 0.13 \; rad, & \mbox{for }27 \leq t < 33 \\ - 0.2\; rad, & \mbox{for }33 \leq t \leq 40 \end{cases} $$ $$ <br\> <br\>
 * $$\displaystyle
 * $$\displaystyle
 * <p style="text-align:right;">$$\displaystyle
 * }
 * }

The solution should be obtained using Hermite-Simpson (HS) algorithm and Newton-Raphson (NR) method. <br\>

Also find $$ J = \int\limits_{t = 0}^ {y(t)dt} $$ where $$ \displaystyle y({t_f}) = 0 $$ for $$ \displaystyle \alpha_1(t) $$ and $$ \displaystyle \alpha_2(t) $$.

Solution
The following code solves the equations of motion listed above for the given control profiles using the Hermite-Simpson with Newton-Raphson algorithm. It contains functions for the main file, dynamics, derivatives, function evaluation, and control profiles.

Figures 7.2.1 and 7.2.2 plot the control inputs to the equations of motion. The left plot in Figure 7.2.1 corresponds to the angle of attack profile specified in the s10 assignment, while the right plot specifies the angle of attack profile for s11. The same thrust profile was used in both cases.
 * HW7.2_Alpha1_Alpha2.png


 * HW7.2_Thrust.png

Figures 7.2.3 - 7.2.6 plot the solutions to the state trajectories found using the Hermite-Simpson and Newton-Raphson algorithm. Again, the left plots correspond to the input profiles from s10 and the right plots correspond to s11. The required integral J from the problem statement corresponds to the area under the curve of the plots in Figure 7.2.6. Figure 7.2.7 plots the horizontal distance versus the vertical distance that the aircraft has traveled (in other words, the spacial representation of the aircraft maneuver).
 * HW7_2(Gamma_t).png


 * HW7_2(V_t).png


 * HW7_2(X_t).png


 * HW7_2(Y_t).png


 * HW7_2(Y_X).png

The integral of y(t) was found for both cases using the code below. The integral for the first angle of attack profile, J1, was calculated to be 16860 m-s, while J2 was found to be 10825 m-s. Since the objective of this problem is to minimize the integral of y(t) over time, the second angle of attack profile (s11) is closer to the optimal solution. The Hermite-Simpson algorithm required 219+1 node points to reach the desired accuracy tolerances. The integrals were found using the the complex trapezoidal rule with values at all of the nodes used.

= Problem 7.3: Solving the Logistic Equation Using Inconsistent Trap-Simpson's rule & Newton-Raphson =

Given: Verhulst Function and Hermite-Simpson Analytical Solution
Refer to lecture slide 41-1 for problem statement

Verhulst's logistic function for population dynamics is given in S10 lecture slide 38-3,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

\dot{x}(t) = f(x) = rx \left(1-\frac{x}{x_{max}}\right).

$$     (7.3.1)
 * <p style="text-align:right">
 * }

For S11, $$\displaystyle x_{max}=15;\ r = 1.4,\ $$ and $$\displaystyle t \in [0,20] $$

The analytical solution is given in S10 lecture slide 39-1,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

x(t) = \frac{x_0x_{max}e^{rt}}{x_{max} + x_0(e^{rt}-1)}.

$$     (7.3.2)
 * <p style="text-align:right">
 * }

Find: Numerical Solution to Verhulst Eqn. Using Inconsistent Trap-Simpson's Rule
To solve logistic equation defined by Eq 7.8.1 for two initial conditions, $$\displaystyle x_0=3\ and\ x_0 = 9,\ $$ using Simpson's rule given by (1) in lecture slide 38-3,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

z_{i+1} = z_{i} + \frac{h/2}{3} \left[f_i + 4f_{(i+1/2)} +f_{i+1}\right].

$$     (7.3.4)
 * <p style="text-align:right">
 * }

The Inconsistent Trap rule for finding $$z_{i+1/2}$$ is given by,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

z_{i+1/2} = \frac{1}{2} [z_{i} + z_{i+1}].

$$     (7.3.5)
 * <p style="text-align:right">
 * }

Solution: Newton-Raphson Method, Solution and Error Plots and Matlab Code
Consider Eq 7.3.1 in which $$\displaystyle f_{(i+1/2)}$$ is a function of $$\displaystyle z_{i}$$ and $$\displaystyle z_{i+1}$$, where $$\displaystyle z_i, z_{i+1} $$ goes to values of $$\displaystyle x$$ at $$\displaystyle t = t_i$$ and $$\displaystyle t = t_{i+1}$$.

Because this is an initial value problem, $$\displaystyle z_i$$ is known. So, Eq 7.3.4 becomes a function of $$\displaystyle z_{i+1}$$ given by,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

F(z_{i+1}) = 0

$$     (7.3.6)
 * <p style="text-align:right">
 * }

We find the root $$\displaystyle z_{i+1}$$ of Eq 7.3.6 using Newton-Raphson method in lecture 26-2,


 * {| style="width:100%" border="0"

$$\displaystyle
 * style="width:95%" |
 * style="width:95%" |

z_{i+1}^{(k+1)} = z_{i+1}^{(k)} - \left[\frac{dF(z_{i+1}^{(k)})}{dz}\right]^{-1} F(z_{i+1}^{(k)})

$$     (7.3.6)
 * <p style="text-align:right">
 * }

The Newton-Raphson iteration, starting with an initial guess as $$\displaystyle z_{i+1}^{0} = z_{i}$$, is stopped once the absolute tolerance


 * {| style="width:100%" border="0"

$$\displaystyle \begin{align} AbsTol &= ||z_{i+1}^{(k+1)} - z_{i+1}^{(k)}|| \\ &\le 10^{-6}\\ \end{align} $$     (7.3.7)
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right">
 * }

is satisfied.

Inconsistent Trap-Simpson Plots
Below are the plots of the Inconsistent Trap-Simpson solution of the Verhulst Function and the error, initial population of 3, h=0.125.

Below are the plots of the Inconsistent Trap-Simpson solution of the Verhulst Function and the error, initial population of 9.

The code for the case with initial population of 9 is very similar to the code post above. The only difference is that "x0=9" instead of "x0=3."

Error Plots of Hermite Simpson Method
Below are the plots of the error of the Hermite Simpson solution compared to the analytical solution of the Verhulst Function.

The code for the case with initial population of 9 is very similar to the code post above. The only difference is that "x0=9" instead of "x0=3."

Conclusion
Comparing the graphs of the errors of the two methods (Inconsistent Trap-Simpson and Hermite Simpson) we see that the error in the Hermite Simpson method (order $$10^{-6}$$) is orders of magnitude less than that in the Inconsistent Trap-Simpson method (order $$10^{-3}$$). This is to be expected because the Hermite Simpson method is consistent in it's picking of the $$ z_{i+\frac{1}{2}} $$ term, whereas the Inconsistent Trap-Simpson is, as its name says, inconsistent.

=References=

=Contributing Members & Referenced Lecture=