This post is established based on the project I have read on the book “A First Course in Mathematical Modeling” by Frank R. Giordano et al. The project is about the spread of a contagious disease. The spread of a contagious disease can be modeled by a system of differential equations.

Problem Background

Consider the following ordinary differential equation model for the spread of a communicable disease:

where $N$ is measured in $100$’s. Analyze the behavior of this differential equation as follows.

1. Since this is an autonomous differential equation, perform a qualitative graphical analysis as discussed in Section 11.4.
   • Plot $dN/dt$ versus $N$. Find and label all rest points (equilibrium points).
   • Estimate the value where the rate of change of the disease is the fastest. Justify your answer.
   • Plot $N$ versus $t$ for each of the following initial conditions: \(N(0) = 2, \quad N(0) = 7, \text{ and } N(0) = 14.\)
   • Describe the stability of each rest point (equilibrium point) in the context of the spread of the disease.
2. Obtain a slope field plot of this differential equation. Briefly analyze the slope field plot. Compare it to your qualitative plot in part 1.
3. Solve this differential equation using separation of variables. Compare your solution to the slope field plot and the qualitative plot in part 1.
4. Compare the time $t$ when $N$ is changing the fastest using the initial condition $N(0) = 2$. Compare your answer to your qualitative plot in part 1.
5. Use Euler’s Method with step sizes of $h=1$ and then $h=0.1$ to approximate the solution to the differential equation for $N(0.5)$ and $N(5)$. Find the relative error at these values. Obtain graphical plots of the numerical solution. How do they compare to the plots of the actual solution?
6. Repeat part 5 using either the improved Euler’s Method or the Runge-Kutta Method.

Preliminary

Before we proceed to the solution, we need to review the key concepts and methods that will be used in solving the problem.

Differential Equations

A differential equation is an equation that relates a function to its derivatives. Differential equations are used to model a wide range of physical phenomena, from population growth to the spread of diseases. There are two main types of differential equations: ordinary differential equations (ODEs) and partial differential equations (PDEs). ODEs involve functions of a single variable, while PDEs involve functions of multiple variables.

Autonomous Differential Equations

An autonomous differential equation is a differential equation where the independent variable does not appear explicitly. In other words, the rate of change of the dependent variable depends only on the value of the dependent variable itself. Autonomous differential equations are often used to model systems that evolve over time without external influences.

Euler’s Method

Euler’s method is a simple numerical method for solving ordinary differential equations. It is based on the idea of approximating the solution by taking small steps along the direction field of the differential equation. Euler’s method is easy to implement and provides a good approximation of the solution for small step sizes.

It is given by the formula:

where $y_n$ is the approximate solution at time $t_n$, $h$ is the step size, and $f(y_n, t_n)$ is the derivative of the function at $(y_n, t_n)$. The method is iterative, with each step updating the solution based on the previous value.

Runge-Kutta Method

The Runge-Kutta method is a family of numerical methods for solving ordinary differential equations. The most common form is the fourth-order Runge-Kutta method (RK4), which provides a good balance between accuracy and computational efficiency. The RK4 method is based on a weighted average of four estimates of the derivative at different points in the interval.

The RK4 method is given by the formula:

where

The RK4 method is more accurate than Euler’s method and provides a good approximation of the solution with fewer steps. We will use Euler’s Method and Runge-Kutta Method to approximate the solution to the differential equation in this problem, and hope you can see how they work.

Slope Field Plot

A slope field plot is a graphical representation of a first-order differential equation. It shows the direction of the solution curves at different points in the $xy$-plane. The slope field is constructed by drawing short line segments at each point, with the slope of the line segment representing the value of the derivative at that point. The slope field provides a visual representation of the behavior of the differential equation and helps in understanding the solutions.

Equilibrium Points

Equilibrium points are points where the derivative of the function is zero. In the context of a differential equation, equilibrium points are points where the rate of change of the dependent variable is zero. Equilibrium points are important in analyzing the behavior of the system and determining the stability of the solutions.

Solution

This is a very typical problem in mathematical modeling. The spread of a contagious disease can be modeled by a system of differential equations. Usually, one may use the SIR model to model the spread of a contagious disease. The SIR model is a simple mathematical model of infectious disease spread. The model divides the population into three compartments: S for the number of susceptible, I for the number of infectious, and R for the number of recovered or immune individuals. But in here we are given a simple differential equation model for the spread of a communicable disease.

Therefore, one may simply use Python to solve this problem. The following code is the implementation of the solution to the problem.

Question 1

First we import the necessary libraries and define the differential equation.

# Importing the required libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve, root

# Setting the font properties
import matplotlib as mpl
mpl.rcParams['mathtext.fontset'] = 'stix'
mpl.rcParams['mathtext.rm'] = 'serif'
mpl.rcParams['font.size'] = 14
mpl.rcParams['font.family'] = 'STIXGeneral'

Then we define the differential equation.

# Define the differential equation
def dNdt(N, t):
    return 0.25 * N * (10 - N)

After that, we have to create arrays for plotting the graphs.

# Create arrays for plotting
N = np.linspace(0, 15, 1000)
dN = dNdt(N, 0)

Now we can plot the graph of $dN/dt$ versus $N$.

# Plot dN/dt vs N
plt.figure(figsize=(10, 6))
plt.plot(N, dN)
plt.axhline(y=0, color='r', linestyle='--')
plt.axvline(x=0, color='g', linestyle='--')
plt.axvline(x=10, color='g', linestyle='--')
plt.title('dN/dt vs N')
plt.xlabel('N')
plt.ylabel('dN/dt')
plt.grid(True)
plt.show()

The graph of $dN/dt$ versus $N$ is shown below.

Then we find the point of fastest rate of change of the disease.

# Find the point of fastest change
N_fastest = N[np.argmax(dN)] # argmax returns the index of the maximum value
print(f"The rate of change is fastest when N = {N_fastest:.5f}")

The output then is

The rate of change is fastest when N = 5.00000

Then solve the differential equation using odeint function.

# Solve the ODE for different initial conditions
t = np.linspace(0, 20, 1000)
initial_conditions = [2, 7, 14]
colors = ['b', 'g', 'r']

plt.figure(figsize=(10, 6))
for N0, color in zip(initial_conditions, colors):
    solution = odeint(dNdt, N0, t)
    plt.plot(t, solution, color, label=f'N(0) = {N0}')

plt.axhline(y=10, color='k', linestyle='--')
plt.title('N vs t for different initial conditions')
plt.xlabel('t')
plt.ylabel('N')
plt.legend()
plt.grid(True)
plt.show()

Then we obtain the following figure.

Time to move on to the next question.

Question 2

Now we have to obtain a slope field plot of this differential equation.

# Create a grid of points
N = np.linspace(0, 15, 20)
t = np.linspace(0, 10, 20)
N, t = np.meshgrid(N, t)

which gives us the grid of points.

Then we calculate the slope at each point.

# Calculate dN/dt at each point
dN = dNdt(N, t)
dt = np.ones_like(N)

Then we normalize the slopes for better visualization.

# Normalize the slopes
length = np.sqrt(dt**2 + dN**2)
dN /= length
dt /= length

Finally, we plot the slope field.

# Create the plot
plt.figure(figsize=(12, 8))
plt.quiver(t, N, dt, dN, pivot='mid', color='b', alpha=0.8)
# Add some solution curves
for N0 in [2, 7, 14]:
    solution = odeint(dNdt, N0, np.linspace(0, 10, 1000))
    plt.plot(np.linspace(0, 10, 1000), solution, 'r-')

plt.title('Slope Field and Solution Curves')
plt.xlabel('t')
plt.ylabel('N')
plt.xlim(0, 10)
plt.ylim(0, 15)
plt.grid(True)
plt.show()

The slope field plot is shown below.

The slope field plot is a visual representation of the differential equation. The slope field plot shows the direction of the solution curves at each point in the $N-t$ plane. The solution curves are the curves that satisfy the differential equation. The slope field plot is a useful tool for understanding the behavior of the differential equation.

Analysis of the slope field plot:

• The arrows in the plot represent the direction and magnitude of dN/dt at each point.
• We can see two horizontal lines where the arrows become very small or disappear completely. These correspond to our equilibrium points at $N = 0$ and $N = 10$.
• Below $N = 10$, the arrows point upward, indicating that $N$ is increasing. Above $N = 10$, the arrows point downward, showing that $N$ is decreasing.
• The arrows are longest around $N = 5$, confirming our earlier finding that the rate of change is fastest near this value.
• Near $N = 0$, the arrows are very small, showing that change is very slow when $N$ is close to zero.

Comparison to the qualitative plot in part 1:

• This slope field plot confirms our analysis from part 1. We can clearly see the stable equilibrium at $N = 10$ and the unstable equilibrium at $N = 0$.
• The solution curves (in red) match our qualitative sketches from part 1. Solutions starting below 10 increase towards $10$, while solutions starting above $10$ decrease towards $10$.
• The slope field provides a more comprehensive view of the behavior of solutions for all possible initial conditions, not just the three we considered in part 1.
• We can see that the rate of change is indeed fastest around $N = 5$, as the arrows are longest in this region.

Question 3

We solve the differential equation using separation of variables. Since it is trivial, we skip this part, one may refer to any ODE textbook for the solution.

Now we plot the solution over the interval of time from $0$ to $10$.

First, we define the analytical solution.

def N_analytical(t):
    return 10 * np.exp(2.5*t) / (4 + np.exp(2.5*t))

Then we redefine the time array

# Redefine the time array
t = np.linspace(0, 10, 1000)

Finally, we plot the analytical solution with comparison to the numerical solution.

# Plot analytical solution
plt.figure(figsize=(12, 8))
plt.plot(t, N_analytical(t), 'b-', label='Analytical Solution')

# Plot numerical solutions for comparison
for N0 in [2, 7, 14]:
    solution = odeint(dNdt, N0, t)
    plt.plot(t, solution, '--', label=f'Numerical Solution (N(0)={N0})')

plt.title('Analytical Solution vs Numerical Solutions')
plt.xlabel('t')
plt.ylabel('N')
plt.legend()
plt.grid(True)
plt.ylim(0, 15)
plt.show()

The plot of the analytical solution versus the numerical solutions is shown below.

As always, we do have analysis of the plot.

Analysis and Comparison:

1. The analytical solution matches very closely with the numerical solution for $N(0) = 2$, which is expected as we used this initial condition to solve for the constant of integration.
2. The analytical solution confirms our qualitative analysis from Part 1:
   • It starts at $N(0) = 2$ and increases towards $10$.
   • The rate of increase is highest in the middle part of the curve.
   • The solution approaches but never quite reaches $N = 10$.
3. Comparing to the slope field plot from Part 2:
   • The solution curve follows the direction of the arrows in the slope field.
   • It shows the fastest change (steepest slope) in the middle region where the arrows were longest in the slope field plot.
4. The analytical solution provides a smooth, continuous curve that fits perfectly with our understanding of the differential equation’s behavior.
5. It confirms that $10$ is indeed a stable equilibrium point, as all solutions approach this value as t increases, regardless of the initial condition.

Question 4

We also redefine the functions

# Define the functions 
def N(t):
    return 10 * np.exp(2.5*t) / (4 + np.exp(2.5*t))

def dNdt(t):
    return 25 * np.exp(2.5*t) * (4 + np.exp(2.5*t)) / (4 + np.exp(2.5*t))**2 - \
           25 * np.exp(5*t) / (4 + np.exp(2.5*t))**2

def d2Ndt2(t):
    return 62.5 * np.exp(2.5*t) * (4 - np.exp(2.5*t)) / (4 + np.exp(2.5*t))**3

Then we calculate the time when $N$ is changing the fastest. One may determine it by finding the time when the second derivative of $N$ is zero.

# Find the time of fastest change
t_fastest = root(d2Ndt2, 0.5).x[0]

print(f"The population is changing fastest at t ≈ {t_fastest:.4f}")
print(f"At this time, N ≈ {N(t_fastest):.4f}")

In here, I use root function from scipy.optimize to find the root of the second derivative of $N$. However, one may use fsolve function as well. Anyway, both of them will give the same result and the initial guess is quite important for both functions.

The output then is

The population is changing fastest at t ≈ 0.5545
At this time, N ≈ 5.0000

The result is consistent with our earlier analysis that the rate of change is fastest when $N = 5$.

We also plot the graph of $N$ versus $t$ for the initial condition $N(0) = 2$.

# Plot N(t) and dN/dt
t = np.linspace(0, 10, 1000)
plt.figure(figsize=(12, 8))
plt.plot(t, N(t), 'b-', label='N(t)')
plt.plot(t, dNdt(t), 'r-', label='dN/dt')
plt.axvline(x=t_fastest, color='g', linestyle='--', label='Fastest change')
plt.axhline(y=N(t_fastest), color='g', linestyle='--')
plt.title('N(t) and dN/dt vs t')
plt.xlabel('t')
plt.ylabel('N / dN/dt')
plt.legend()
plt.grid(True)
plt.show()

The plot of $N(t)$ and $dN/dt$ versus $t$ is shown below.

Analysis and Comparison to Part 1:

• The population is changing fastest at $t ≈ 0.5545$, when $N ≈ 5.0000$.
• This confirms our qualitative analysis from Part 1, where we estimated that the rate of change would be fastest when $N = 5$.
• The plot shows that $dN/dt$ (the red curve) reaches its maximum at the same point where $N(t)$ (the blue curve) has its steepest slope.
• This point occurs early in the growth process, when the population has reached about half of its carrying capacity $(10)$.
• After this point, the rate of change $(dN/dt)$ starts to decrease, which is reflected in the decreasing slope of $N(t)$.
• This behavior aligns with our slope field plot from Part 2, where we observed that the arrows (representing $dN/dt$) were longest around $N = 5$.

Question 5

We use Euler’s Method to approximate the solution to the differential equation for $N(0.5)$ and $N(5)$. We start by defining the Euler’s Method function.

# Define the functions
def dNdt(N, t):
    return 0.25 * N * (10 - N)

def exact_solution(t):
    return 10 * np.exp(2.5*t) / (4 + np.exp(2.5*t))

# Solve the ODE using Euler's method
def euler_method(f, y0, t0, tf, h):
    t = np.arange(t0, tf+h, h)
    y = np.zeros(len(t))
    y[0] = y0
    for i in range(1, len(t)):
        y[i] = y[i-1] + h * f(y[i-1], t[i-1])
    return t, y

Then we set up the parameters for the Euler’s Method.

# Parameters
N0 = 2
t0 = 0
tf = 5

Now we approximate the solution using Euler’s Method with step sizes of $h=1$ and $h=0.1$.

# Euler's method with h=1
t1, N1 = euler_method(dNdt, N0, t0, tf, 1)

# Euler's method with h=0.1
t01, N01 = euler_method(dNdt, N0, t0, tf, 0.1)

We also calculate the exact solution for comparison.

# Exact solution
t_exact = np.linspace(t0, tf, 1000)
N_exact = exact_solution(t_exact)

Moreover, we calculate the relative error at $N(0.5)$ and $N(5)$.

# Calculate approximations and errors
N05_h1 = np.interp(0.5, t1, N1)
N5_h1 = N1[-1]
N05_h01 = np.interp(0.5, t01, N01)
N5_h01 = N01[-1]

N05_exact = exact_solution(0.5)
N5_exact = exact_solution(5)

rel_error_05_h1 = abs((N05_h1 - N05_exact) / N05_exact) * 100
rel_error_5_h1 = abs((N5_h1 - N5_exact) / N5_exact) * 100
rel_error_05_h01 = abs((N05_h01 - N05_exact) / N05_exact) * 100
rel_error_5_h01 = abs((N5_h01 - N5_exact) / N5_exact) * 100

print(f"N(0.5) approximation (h=1): {N05_h1:.4f}, Relative Error: {rel_error_05_h1:.2f}%")
print(f"N(5) approximation (h=1): {N5_h1:.4f}, Relative Error: {rel_error_5_h1:.2f}%")
print(f"N(0.5) approximation (h=0.1): {N05_h01:.4f}, Relative Error: {rel_error_05_h01:.2f}%")
print(f"N(5) approximation (h=0.1): {N5_h01:.4f}, Relative Error: {rel_error_5_h01:.2f}%")

The output then is

The rate of change is fastest when N = 5.00000
The population is changing fastest at t ≈ 0.5545
At this time, N ≈ 5.0000

N(0.5) approximation (h=1): 4.0000, Relative Error: 14.16%
N(5) approximation (h=1): 6.0000, Relative Error: 40.00%
N(0.5) approximation (h=0.1): 4.5204, Relative Error: 2.99%
N(5) approximation (h=0.1): 10.0000, Relative Error: 0.00%

Finally, we plot the numerical solution obtained using Euler’s Method.

# Plotting
plt.figure(figsize=(12, 8))
plt.plot(t_exact, N_exact, 'b-', label='Exact Solution')
plt.plot(t1, N1, 'ro-', label="Euler's Method (h=1)")
plt.plot(t01, N01, 'g*-', label="Euler's Method (h=0.1)")
plt.title("Comparison of Euler's Method with Exact Solution", fontsize=18)
plt.xlabel('t')
plt.ylabel('N')
plt.legend()
plt.grid(True)
plt.show()

The plot of the numerical solution obtained using Euler’s Method is shown below.

Analysis of Results:

1. Approximations and Relative Errors:
   • For $h=1$: $N(0.5) ≈ 4.0000$, Relative Error: $14.16%$, $N(5) ≈ 6.0000$, Relative Error: $40.00%$
   • For $h=0.1$: $N(0.5) ≈ 4.5204$, Relative Error: $2.99%$, $N(5) ≈ 10.00$, Relative Error: $0.00%$
2. Observations:
   • The smaller step size ($h=0.1$) provides much more accurate approximations than the larger step size ($h=1$) which is expected when using Euler’s method.
   • The relative errors for $h=0.1$ are significantly smaller than for $h=1$.
   • The approximations for $N(5)$ are generally more accurate than for $N(0.5)$, likely because the solution is approaching the stable equilibrium point at $N=10$.
3. Graphical Comparison:
   • The plot shows that Euler’s method with $h=1$ (red line) provides a rough approximation of the exact solution (blue line).
   • Euler’s method with $h=0.1$ (green line) provides a much closer approximation to the exact solution.
   • The $h=1$ approximation underestimates the solution for most of the interval, especially in the early stages of growth.
   • The $h=0.1$ approximation is nearly indistinguishable from the exact solution on this scale.
4. Limitations of Euler’s Method:
   • Euler’s method tends to underestimate the solution for this problem because it uses a linear approximation for each step.
   • The larger step size ($h=1$) leads to significant accumulation of error over time.
   • Even with $h=0.1$, there is still some error, although it’s much smaller.
5. Comparison to earlier parts:
   • These numerical results confirm our earlier analysis about the general behavior of the solution.
   • The solution starts at $N(0)=2$ and approaches but doesn’t quite reach $N=10$ by $t=5$.
   • The steepest growth occurs in the early-to-middle part of the time interval, as we observed in earlier parts.

Question 6

Finally, we repeat the process by using Runge-Kutta Method. We start by defining the Runge-Kutta Method function.

# Implementing the Runge-Kutta 4th order method
def runge_kutta_4(f, y0, t0, tf, h):
    t = np.arange(t0, tf+h, h)
    y = np.zeros(len(t))
    y[0] = y0
    for i in range(1, len(t)):
        k1 = h * f(y[i-1], t[i-1])
        k2 = h * f(y[i-1] + 0.5*k1, t[i-1] + 0.5*h)
        k3 = h * f(y[i-1] + 0.5*k2, t[i-1] + 0.5*h)
        k4 = h * f(y[i-1] + k3, t[i-1] + h)
        y[i] = y[i-1] + (k1 + 2*k2 + 2*k3 + k4) / 6
    return t, y

Then we also set up the parameters for the Runge-Kutta Method.

# Parameters
N0 = 2
t0 = 0
tf = 5

Now we estimate the solution using Runge-Kutta Method with step sizes of $h=1$ and $h=0.1$.

# RK4 with h=1
t1, N1 = runge_kutta_4(dNdt, N0, t0, tf, 1)

# RK4 with h=0.1
t01, N01 = runge_kutta_4(dNdt, N0, t0, tf, 0.1)

# Exact solution
t_exact = np.linspace(t0, tf, 1000)
N_exact = exact_solution(t_exact)

We also calculate the relative error at $N(0.5)$ and $N(5)$.

# Calculate approximations and errors
N05_h1 = np.interp(0.5, t1, N1)
N5_h1 = N1[-1]
N05_h01 = np.interp(0.5, t01, N01)
N5_h01 = N01[-1]

N05_exact = exact_solution(0.5)
N5_exact = exact_solution(5)

rel_error_05_h1 = abs((N05_h1 - N05_exact) / N05_exact) * 100
rel_error_5_h1 = abs((N5_h1 - N5_exact) / N5_exact) * 100
rel_error_05_h01 = abs((N05_h01 - N05_exact) / N05_exact) * 100
rel_error_5_h01 = abs((N5_h01 - N5_exact) / N5_exact) * 100

print(f"N(0.5) approximation (h=1): {N05_h1:.4f}, Relative Error: {rel_error_05_h1:.4f}%")
print(f"N(5) approximation (h=1): {N5_h1:.4f}, Relative Error: {rel_error_5_h1:.4f}%")
print(f"N(0.5) approximation (h=0.1): {N05_h01:.4f}, Relative Error: {rel_error_05_h01:.4f}%")
print(f"N(5) approximation (h=0.1): {N5_h01:.4f}, Relative Error: {rel_error_5_h01:.4f}%")

The output then is

# Calculate approximations and errors
N05_h1 = np.interp(0.5, t1, N1)
N5_h1 = N1[-1]
N05_h01 = np.interp(0.5, t01, N01)
N5_h01 = N01[-1]

N05_exact = exact_solution(0.5)
N5_exact = exact_solution(5)

rel_error_05_h1 = abs((N05_h1 - N05_exact) / N05_exact) * 100
rel_error_5_h1 = abs((N5_h1 - N5_exact) / N5_exact) * 100
rel_error_05_h01 = abs((N05_h01 - N05_exact) / N05_exact) * 100
rel_error_5_h01 = abs((N5_h01 - N5_exact) / N5_exact) * 100

print(f"N(0.5) approximation (h=1): {N05_h1:.4f}, Relative Error: {rel_error_05_h1:.4f}%")
print(f"N(5) approximation (h=1): {N5_h1:.4f}, Relative Error: {rel_error_5_h1:.4f}%")
print(f"N(0.5) approximation (h=0.1): {N05_h01:.4f}, Relative Error: {rel_error_05_h01:.4f}%")
print(f"N(5) approximation (h=0.1): {N5_h01:.4f}, Relative Error: {rel_error_5_h01:.4f}%")

Finally, we plot the numerical solution obtained using Runge-Kutta Method.

# Plotting# Plotting
plt.figure(figsize=(12, 8))
plt.plot(t_exact, N_exact, 'b-', label='Exact Solution')
plt.plot(t1, N1, 'ro-', label="RK4 (h=1)")
plt.plot(t01, N01, 'g.-', label="RK4 (h=0.1)")
plt.title("Comparison of RK4 Method with Exact Solution")
plt.xlabel('t')
plt.ylabel('N')
plt.legend()
plt.grid(True)
plt.show()

The plot of the numerical solution obtained using Runge-Kutta Method is shown below.

Analysis of Results:

1. Approximations and Relative Errors:
   • For $h=1$: $N(0.5) ≈ 4.6758$, Relative Error: $0.3432%$; $N(5) ≈ 9.7076$, Relative Error: $2.9225%$
   • For $h=0.1$: $N(0.5) ≈ 4.6598$, Relative Error: $0.0003%$; $N(5) ≈ 9.9999$, Relative Error: $0.0000%$
2. Comparison with Euler’s Method:
   • RK4 is significantly more accurate than Euler’s method, even with the larger step size.
   • The relative errors for RK4 are much smaller than those for Euler’s method, especially for $h=1$.
3. Observations:
   • The RK4 method with $h=1$ provides very accurate results, with relative errors less than $0.5%$.
   • The RK4 method with $h=0.1$ provides extremely accurate results, with relative errors less than $0.0001%$.
   • Both step sizes give excellent approximations, but $h=0.1$ is notably more accurate.
4. Graphical Comparison:
   • The plot shows that both RK4 approximations ($h=1$ and $h=0.1$) are visually indistinguishable from the exact solution.
   • This is a significant improvement over Euler’s method, where the $h=1$ approximation was noticeably different from the exact solution.
5. Advantages of RK4:
   • RK4 provides much higher accuracy than Euler’s method, even with larger step sizes.
   • This allows for efficient computation with fewer steps while maintaining high accuracy.
   • The method captures the curvature of the solution much better than Euler’s method, which only uses linear approximations.
6. Comparison to earlier parts:
   • The RK4 results confirm our earlier analyses about the behavior of the solution.
   • The high accuracy of RK4 gives us confidence in our understanding of the solution’s behavior throughout the time interval.

Conclusion

In this post, we have analyzed the spread of a contagious disease using a simple differential equation model. We performed a qualitative graphical analysis, obtained a slope field plot, solved the differential equation using separation of variables, and compared the results to the slope field plot. We also used Euler’s method and Runge-Kutta method to approximate the solution to the differential equation and compared the results to the exact solution. The results showed that the Runge-Kutta method provided much higher accuracy than Euler’s method, even with larger step sizes. The RK4 method was able to capture the curvature of the solution much better than Euler’s method, leading to more accurate approximations. The numerical results confirmed our earlier analyses about the behavior of the solution and provided valuable insights into the spread of the contagious disease.