Applied Mathematics for Engineering

4. Numerical Methods

Numerical methods are essential tools in engineering for solving mathematical problems when analytical solutions are not feasible or computationally efficient. In this article, we will discuss some common numerical methods used in engineering, including numerical integration, numerical differentiation, and solving ordinary and partial differential equations numerically. We will also demonstrate how to implement these methods in Python using the SciPy library.

4.1. Numerical Integration

Numerical integration, also known as quadrature, is the process of computing the definite integral of a function. It is often used to find the area under a curve or the value of a definite integral when an analytical solution is challenging to obtain.

One popular method for numerical integration is the Simpson’s rule. The rule approximates the definite integral of a function by dividing the integration interval into subintervals and approximating the function with quadratic polynomials over each subinterval. The formula for Simpson’s rule is:

\[\int_a^b f(x) dx \approx \frac{b-a}{6n} [f(a) + 4\sum_{i=1}^{n} f(x_{2i-1}) + 2\sum_{i=1}^{n-1} f(x_{2i}) + f(b)]\]

In mechatronics engineering, numerical integration is used in various applications, such as computing the work done by a force on a robot arm, estimating the power consumed by an electric motor, or calculating the response of a control system.

4.1.1. Using SciPy for Numerical Integration

To perform numerical integration in Python, you can use the scipy.integrate.quad function:

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.sin(x)

result, error = quad(func, 0, np.pi)
print(f"Result: {result}, Error: {error}")

4.2. Numerical Differentiation

Numerical differentiation is the process of estimating the derivative of a function using discrete data points. Finite difference methods are common techniques for numerical differentiation, with forward, backward, and central difference methods being the most popular.

The central difference method is more accurate than the forward and backward difference methods. The formula for the central difference method is:

\[f'(x_i) \approx \frac{f(x_{i+1}) - f(x_{i-1})}{2h}\]

where \(h\) is the step size.

In engineering, numerical differentiation is applied to estimate the rate of change of various quantities, such as the velocity and acceleration of a robotic system or the rate of change of an electrical current in a circuit.

4.2.1. Using SciPy for Numerical Differentiation

To perform numerical differentiation in Python, you can use the scipy.misc.derivative function:

import numpy as np
from scipy.misc import derivative

def func(x):
    return np.sin(x)

result = derivative(func, np.pi/4, dx=1e-6)
print(f"Result: {result}")

4.3. Solving ODEs and PDEs Numerically

Ordinary differential equations (ODEs) and partial differential equations (PDEs) are common in engineering problems. While some ODEs and PDEs have analytical solutions, many require numerical methods for their solutions.

The Runge-Kutta method is a widely-used numerical method for solving ODEs. The fourth-order Runge-Kutta method (RK4) is particularly popular due to its accuracy and simplicity.

In engineering, solving ODEs and PDEs numerically is crucial for designing and analyzing systems such as control systems, robotic systems, and dynamic systems. For example, you may need to solve ODEs to simulate the trajectory of a mobile robot or the response of a feedback control system, while PDEs are commonly used to model heat transfer or fluid dynamics in various applications.

4.3.1. Using SciPy for Solving ODEs with the Runge-Kutta Method

To solve ODEs in Python, you can use the scipy.integrate.solve_ivp function. The following example demonstrates how to solve a simple first-order ODE using the Runge-Kutta method:

import numpy as np
from scipy.integrate import solve_ivp

def func(t, y):
    return np.cos(t)

solution = solve_ivp(func, [0, 10], [0], method='RK45', t_eval=np.linspace(0, 10, 100))

import matplotlib.pyplot as plt

plt.plot(solution.t, solution.y[0], label='Numerical solution')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()

This example solves the ODE \(y'(t) = \cos(t)\) with the initial condition \(y(0) = 0\) using the Runge-Kutta method.

4.4. Exercises

Example 1

Compute the numerical integral of the function \(f(x) = e^{-x^2}\) on the interval \(0 \le x \le 3\) using the Simpson’s rule. Provide the result and the error.

Solution:

import numpy as np
from scipy.integrate import quad

def func(x):
    return np.exp(-x**2)

result, error = quad(func, 0, 3)
print(f"Result: {result}, Error: {error}")

Example 2

Compute the numerical derivative of the function \(f(x) = e^x\) at the point \(x = 2\) using the central difference method. Provide the result.

Solution:

import numpy as np
from scipy.misc import derivative

def func(x):
    return np.exp(x)

result = derivative(func, 2, dx=1e-6)
print(f"Result: {result}")

Example 3

Solve the following ODE using the Runge-Kutta method:

\(y''(t) + y'(t) + y(t) = \sin(t)\), with initial conditions \(y(0) = 1\) and \(y'(0) = 0\).

Solution:

import numpy as np
from scipy.integrate import solve_ivp

def func(t, y):
    return [y[1], np.sin(t) - y[1] - y[0]]

solution = solve_ivp(func, [0, 10], [1, 0], method='RK45', t_eval=np.linspace(0, 10, 100))

import matplotlib.pyplot as plt

plt.plot(solution.t, solution.y[0], label='Numerical solution')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()

Example 4

How can the Runge-Kutta method be applied to simulate the motion of a simple pendulum in Python, and what are its advantages over other numerical methods?

Solution:

The simple pendulum is a classic example of a physical system that can be modeled using ODEs. The Runge-Kutta method is a powerful numerical integration technique that can be used to simulate the motion of the pendulum in Python. Here is an example of implementing the 4th order Runge-Kutta method to simulate the motion of a simple pendulum.

The motion of a simple pendulum can be described by the following second-order ODE:

\[\frac{d^2 \theta}{dt^2} + \frac{g}{l} \sin{\theta} = 0\]

where \(\theta\) is the angle of the pendulum, \(g\) is the gravitational constant, and \(l\) is the length of the pendulum. We can convert this second-order ODE into two first-order ODEs by introducing a new variable \(\omega = \frac{d\theta}{dt}\):

\[\begin{align*} \frac{d \theta}{dt} &= \omega \\ \frac{d \omega}{dt} &= -\frac{g}{l} \sin{\theta} \end{align*}\]

We can then implement the 4th order Runge-Kutta method in Python to numerically integrate these ODEs and simulate the motion of the pendulum. Here is an example program:

import numpy as np
import matplotlib.pyplot as plt

# Define the ODEs
def pendulum_ode(t, y, g, l):
    theta, omega = y
    dtheta_dt = omega
    domega_dt = -(g / l) * np.sin(theta)
    return [dtheta_dt, domega_dt]

# Define the simulation parameters
g = 9.81
l = 1.0
t = np.linspace(0, 10, 1000)
y0 = [np.pi/4, 0]

# Use the Runge-Kutta method to simulate the motion of the pendulum
sol = np.zeros((len(t), 2))
sol[0] = y0
for i in range(len(t) - 1):
    h = t[i+1] - t[i]
    k1 = h * np.array(pendulum_ode(t[i], sol[i], g, l))
    k2 = h * np.array(pendulum_ode(t[i] + h/2, sol[i] + k1/2, g, l))
    k3 = h * np.array(pendulum_ode(t[i] + h/2, sol[i] + k2/2, g, l))
    k4 = h * np.array(pendulum_ode(t[i] + h, sol[i] + k3, g, l))
    sol[i+1] = sol[i] + (k1 + 2*k2 + 2*k3 + k4) / 6

# Plot the results
plt.plot(t, sol[:, 0])
plt.xlabel('Time (s)')
plt.ylabel('Angle (rad)')
plt.show()

In this program, we define the ODEs representing the motion of the simple pendulum, and we implement the 4th order Runge-Kutta method to simulate the motion over time. We then plot the angle of the pendulum as a function of time.

The Runge-Kutta method provides higher order approximations and allows for smaller time steps, resulting in more accurate solutions. It is also more robust and stable for stiff ODEs compared to other numerical methods like the Euler method or the Midpoint method. These advantages make the Runge-Kutta method a popular choice for simulating complex systems in science and engineering.

Example 5

Write a Python program to solve the following one-dimensional heat equation using the finite difference method:

\(\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\),

where \(u(x,t)\) represents the temperature distribution in a rod of length \(L\) at time \(t\). Use the Crank-Nicolson method for time-stepping. The boundary conditions and initial conditions are given by:

\(u(0, t) = u(L, t) = 0\) (Dirichlet boundary conditions), and \(u(x, 0) = f(x)\) for \(0 \le x \le L\),

where \(f(x)\) is an initial temperature distribution function. Consider the following parameters: \(\alpha\) is the diffusion coefficient, \(L\) is the length of the rod, and \(T\) is the final time. Use \(0 \le t \le T\) for the time interval.

Solution:

Here’s an example of Python program that solves the one-dimensional heat equation using the finite difference method with the Crank-Nicolson scheme.

In this program, the solve_heat_equation function takes the parameters L (length of the rod), T (final time), alpha (diffusion coefficient), N (number of spatial grid points), M (number of time steps), and the initial temperature distribution function f(x). It returns the spatial grid points x and the temperature distribution matrix U.

The program solves the heat equation by discretizing the domain, constructing the coefficient matrices A and B, initializing the solution matrix U, and iteratively solving the system of equations using the numpy linalg.solve function. Finally, it plots the temperature distribution at various time levels.

You can modify the example usage and the initial temperature distribution function to suit your specific problem. Please feel free to edit or make contribution to this code on GitHub.

 1import numpy as np
 2import matplotlib.pyplot as plt
 3
 4def solve_heat_equation(L, T, alpha, N, M, f):
 5    # Parameters
 6    h = L / N
 7    k = T / M
 8
 9    # Create grid points
10    x = np.linspace(0, L, N+1)
11
12    # Create initial temperature distribution
13    u0 = f(x)
14
15    # Create coefficient matrices
16    A = np.zeros((N-1, N-1))
17    B = np.zeros((N-1, N-1))
18
19    # Populate the coefficient matrices
20    A[0, 0] = 1 + alpha * k / (2 * h**2)
21    A[0, 1] = -alpha * k / (2 * h**2)
22    B[0, 0] = 1 - alpha * k / (2 * h**2)
23    B[0, 1] = alpha * k / (2 * h**2)
24
25    for i in range(1, N-2):
26        A[i, i-1] = -alpha * k / (2 * h**2)
27        A[i, i] = 1 + alpha * k / h**2
28        A[i, i+1] = -alpha * k / (2 * h**2)
29
30        B[i, i-1] = alpha * k / (2 * h**2)
31        B[i, i] = 1 - alpha * k / h**2
32        B[i, i+1] = alpha * k / (2 * h**2)
33
34    A[N-2, N-3] = -alpha * k / (2 * h**2)
35    A[N-2, N-2] = 1 + alpha * k / h**2
36    B[N-2, N-3] = alpha * k / (2 * h**2)
37    B[N-2, N-2] = 1 - alpha * k / h**2
38
39    # Initialize solution matrix
40    U = np.zeros((N+1, M+1))
41    U[:, 0] = u0
42
43    # Solve the system iteratively
44    for j in range(M):
45        # Use matrix solver to find U at the next time step
46        U[1:N, j+1] = np.linalg.solve(A, np.dot(B, U[1:N, j]))
47
48    return x, U
49
50# Example usage
51L = 1.0
52T = 0.5
53alpha = 0.1
54N = 100  # Number of spatial grid points
55M = 1000  # Number of time steps
56
57# Define initial temperature distribution function
58def f(x):
59    # Heat at one end (e.g., candle)
60    heat_location = 0.2 * L  # Specify the location of heat
61    heat_strength = 10.0  # Adjust the strength of heat
62
63    # Calculate the initial temperature distribution
64    initial_temperature = np.exp(-((x - heat_location) ** 2) / 0.01)
65    initial_temperature *= heat_strength
66
67    return initial_temperature
68
69x, U = solve_heat_equation(L, T, alpha, N, M, f)
70
71# Plot the results using color map and annotations
72plt.figure()
73X, T = np.meshgrid(x, np.linspace(T, 0, M+1))
74plt.pcolormesh(X, T, U.T, shading='auto', cmap='hot')
75plt.colorbar(label='Temperature')
76plt.xlabel('x')
77plt.ylabel('Time')
78plt.title('Temperature Distribution')
79plt.grid(True)
80plt.savefig("1d-heat-distribution.png")
81plt.show()

This example solves the heat equation using the Crank-Nicolson method with given boundary and initial conditions. The solution shows the temperature distribution as a function of position and time.

The f(x) function has been modified to represent the initial temperature distribution where one end of the rod (specified by the heat_location parameter) is exposed to heat, simulating the effect of a candle or any localized heat source. The heat_strength parameter controls the intensity of the heat applied.

The temperature distribution where one end of the rod is hotter than the rest. The color map and annotations illustrate the temperature progression over time, with the heat gradually spreading along the rod.

Feel free to adjust the heat_location and heat_strength parameters to match your desired scenario and explore different initial temperature distributions. You may edit or make contribution to this code on GitHub.

Applied Mathematics for Engineering

4. Numerical Methods

4.1. Numerical Integration

4.1.1. Using SciPy for Numerical Integration

4.2. Numerical Differentiation

4.2.1. Using SciPy for Numerical Differentiation

4.3. Solving ODEs and PDEs Numerically

4.3.1. Using SciPy for Solving ODEs with the Runge-Kutta Method

4.4. Exercises

Add Comment

Comments (1)