Scientific Programming and Visualization

Introduction

Engineering analysis increasingly depends on numerical computation. A mechanical engineer may need to solve simultaneous equations, process sensor measurements, evaluate a model over thousands of operating conditions, or communicate a result through a clear graph. Scientific programming combines programming, numerical methods, data structures, and visual communication to perform these tasks reproducibly.

This chapter uses Python concepts and the NumPy and Matplotlib libraries. The principles also apply to MATLAB, Julia, R, and other scientific-computing environments. The objective is not merely to produce a plot, but to create a defensible computational argument: define the problem, represent the data correctly, compute with appropriate methods, visualize the result honestly, and explain its engineering meaning.

Learning Objectives

By the end of this chapter, students should be able to:

1. Distinguish Python lists from numerical arrays and explain why array operations are useful.
2. Create, index, reshape, and combine vectors and matrices using NumPy.
3. Perform matrix arithmetic and solve linear systems while checking dimensions and units.
4. Generate two-dimensional and multidimensional visualizations with Matplotlib.
5. Design publication-quality graphs with suitable labels, scales, legends, and annotations.
6. Represent uncertainty and communicate trends without misleading the reader.
7. Develop a reproducible engineering analysis and evaluate its limitations.

1. Numerical Data Structures

1.1 Scalars, vectors, and matrices

A scalar is a single quantity, such as temperature or density. A vector is an ordered one-dimensional collection of values. For example, a velocity vector in Cartesian coordinates may be written as

where each component has units of metres per second. A matrix is a rectangular array of values. A matrix may represent a system of equations, a coordinate transformation, or a collection of experimental observations.

In NumPy, the central object is an ndarray:

import numpy as np

x = np.array([0.0, 1.0, 2.0, 3.0])
A = np.array([[2.0, 1.0], [1.0, 3.0]])

The attributes x.shape, x.ndim, and x.dtype describe the array's dimensions, number of dimensions, and data type. A two-dimensional array with shape (m, n) has m rows and n columns.

A Python list stores objects and does not automatically perform mathematical vector operations. For example, [1, 2, 3] 2 repeats the list, whereas np.array([1, 2, 3]) 2 gives [2, 4, 6]. This distinction is important in engineering programs.

1.2 Array construction and indexing

Useful constructors include np.zeros, np.ones, np.arange, and np.linspace. The function linspace(a, b, n) creates n equally spaced values including both endpoints, which is often preferable for plotting a continuous model. Indexing begins at zero:

t = np.linspace(0.0, 10.0, 501)
y = np.sin(t)
first_value = y[0]
last_five = y[-5:]

For a matrix A, A[i, j] selects row i and column j. A[:, 0] selects the first column and A[1, :] selects the second row. A slice usually creates a view rather than an independent copy; use .copy() when later changes must not affect the original array.

1.3 Elementwise and matrix operations

NumPy operators normally act element by element. If a and b have compatible shapes, a * b is elementwise multiplication. Matrix multiplication is written a @ b.

For a matrix A and vector x, the product A @ x is defined when the number of columns of A equals the length of x. If A has shape (m, n), then x must have shape (n,) or (n, 1), and the result has m components.

The dot product of two equal-length vectors is

Forces and displacements illustrate its meaning. If force F is measured in newtons and displacement s in metres, mechanical work is

with units of joules, provided the force is constant over the displacement. In Python, np.dot(F, s) or F @ s can calculate this quantity.

Broadcasting allows NumPy to combine arrays with compatible but unequal shapes. For example, adding a scalar temperature offset to every element of a temperature array is valid. Broadcasting is powerful, but it can hide an unintended shape mismatch; therefore, inspect shapes explicitly during development.

2. Matrix Methods for Engineering Models

2.1 Linear systems

Many engineering models can be expressed as

where A is a coefficient matrix, x is the unknown vector, and b is a known right-hand-side vector. The equation is applicable when the equations are linear in the unknowns. If A is square and nonsingular, a unique solution exists.

Use np.linalg.solve(A, b) rather than explicitly calculating A⁻¹b. Direct inversion is generally less efficient and can amplify numerical errors. For example:

A = np.array([[4.0, 1.0], [2.0, 3.0]])
b = np.array([9.0, 13.0])
x = np.linalg.solve(A, b)
residual = A @ x - b

The solution is x = [1.4, 3.4]. Substitution verifies the result:

4(1.4) + 1(3.4) = 9 and 2(1.4) + 3(3.4) = 13.

The residual vector is

and should be close to zero relative to the scale of b. A small residual does not necessarily guarantee an accurate solution if the matrix is ill-conditioned. The condition number, obtained for example from np.linalg.cond(A), indicates sensitivity to perturbations. A large condition number signals that measurement or rounding errors may be strongly amplified.

2.2 Norms and distances

The Euclidean norm of a vector is

It measures vector magnitude when components use compatible units. In NumPy, np.linalg.norm(x) computes this norm by default. Norms are useful for residual checks, error measures, and geometric distances. Components with different physical units should not be combined without a meaningful nondimensionalization or scaling procedure.

2.3 Numerical differentiation and integration

For tabulated data, a forward-difference approximation to the derivative is

This approximation assumes that ∆x is sufficiently small and that the data are reasonably smooth. Central differences generally provide greater accuracy for interior points:

The trapezoidal approximation to an integral is

It assumes that the function is adequately represented between adjacent samples. In current NumPy versions, np.trapezoid is preferred where available; np.trapz is found in older code. Numerical differentiation tends to amplify measurement noise, whereas integration often smooths it. This difference matters when processing experimental signals.

3. From Model to Computation

A reliable scientific program follows a sequence:

1. State the physical question and define inputs and outputs.
2. Choose a mathematical model and record assumptions.
3. Select units and convert all inputs consistently.
4. Represent the data with suitable array shapes.
5. Compute results using tested numerical operations.
6. Check dimensions, limiting cases, residuals, and conservation laws.
7. Visualize the result and compare it with expected physical behaviour.
8. Save code, input data, software versions, and figures.

For instance, the ideal-gas relation is

where p is absolute pressure in pascals, V is volume in cubic metres, n is amount of substance in moles, R is the gas constant, and T is absolute temperature in kelvins. The model assumes ideal-gas behaviour. Using Celsius directly for T would violate the model's unit requirement and produce an incorrect result.

A program should separate data preparation, calculation, and plotting into functions. This improves testing and allows the same calculation to be used with several datasets.

4. Visualization of Engineering Data

4.1 The grammar of a scientific graph

A graph should answer a question. The horizontal axis normally contains the independent variable and the vertical axis the dependent variable, although this is not universal. Every axis requires a descriptive label and unit. A legend is necessary only when multiple distinguishable series are shown. Titles should identify the subject; explanatory detail can be placed in the caption.

Common plot types include:

A line connecting points implies an ordering or continuous relationship. It should not be used to connect unrelated categories merely because they are listed in sequence.

4.2 Basic Matplotlib workflow

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(6.4, 4.2), constrained_layout=True)
ax.plot(t, y, color='tab:blue', linewidth=2.0, label='Model')
ax.scatter(t_data, y_data, color='black', s=18, label='Measurements')
ax.set_xlabel('Time, t (s)')
ax.set_ylabel('Displacement, x (mm)')
ax.grid(True, alpha=0.3)
ax.legend(frameon=False)
fig.savefig('displacement.png', dpi=300)
plt.show()

The object-oriented interface, using fig and ax, is preferable to relying exclusively on global plotting commands. It scales better to multiple panels and makes the code easier to maintain.

Use dpi=300 or greater for raster output intended for reports, and use PDF or SVG when vector output is suitable. A publication-quality figure also requires readable font sizes, sufficient line contrast, consistent units, and an informative caption. Avoid excessive decoration, three-dimensional effects that obscure values, and colours that cannot be distinguished by readers with colour-vision deficiencies.

4.3 Multidimensional results

Suppose temperature is calculated over a rectangular grid of x and y coordinates. A mesh can be created with np.meshgrid:

x = np.linspace(0, 1, 101)
y = np.linspace(0, 0.5, 51)
X, Y = np.meshgrid(x, y)
T = 20 + 80 * X * (1 - 2 * Y)

Here X, Y, and T have the same shape. A contour plot emphasizes level curves, while a pseudocolour plot emphasizes local magnitude:

fig, ax = plt.subplots(figsize=(6.5, 4.5), constrained_layout=True)
levels = np.linspace(T.min(), T.max(), 12)
cs = ax.contourf(X, Y, T, levels=levels, cmap='viridis')
fig.colorbar(cs, ax=ax, label='Temperature (°C)')
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')

The colour map must be appropriate for the data. Sequential maps suit quantities increasing from low to high. Diverging maps suit data with a meaningful central reference, such as temperature difference from a baseline. A colour bar is essential because colour otherwise has no quantitative interpretation.

For data with a large dynamic range, a logarithmic axis may be appropriate. It should be used only when zero and negative values are not part of the plotted quantity and when multiplicative changes are meaningful. Logarithmic presentation can reveal power-law behaviour but may conceal absolute differences.

5. Uncertainty and Statistical Communication

No measurement is exact. If a measured quantity is reported as x ± u, u should be defined, for example as a standard uncertainty, a confidence interval, or an instrument tolerance. These meanings are not interchangeable.

For independent variables x and y and a calculated quantity q = f(x, y), first-order uncertainty propagation is

This approximation assumes small uncertainties, a sufficiently smooth function, and negligible covariance. With covariance included, the two-variable expression becomes

If x and y are independent, cov(x,y) = 0. For a mean of n independent observations with standard deviation s, the standard error of the mean is

This describes uncertainty in the estimated mean under the independence assumption; it is not the spread of individual observations.

Error bars should state what they represent. They may show instrument uncertainty, standard deviation, standard error, or a confidence interval. A graph with error bars but no definition is scientifically incomplete. When uncertainty is asymmetric, use separate lower and upper errors rather than forcing symmetry.

A shaded band can display model uncertainty across a continuous curve:

ax.plot(t, mean, label='Predicted mean')
ax.fill_between(t, lower, upper, alpha=0.2, label='Uncertainty interval')

Uncertainty also arises from model assumptions, calibration, numerical discretization, and omitted physics. A narrow statistical interval does not imply that the engineering model is accurate if systematic error dominates.

6. A Detailed Applied Case: Thermal Test of an Aluminium Plate

The following is an illustrative case. It demonstrates a complete workflow rather than reporting an actual experiment.

Background and problem

A laboratory team heats one edge of an aluminium plate and records temperature at several positions after 600 seconds. The engineering question is whether the measured temperature distribution is consistent with a simple linear approximation and where additional thermal insulation might be justified.

The positions are x = [0.00, 0.05, 0.10, 0.15, 0.20] m. Measured temperatures are [82, 74, 66, 58, 51] °C, with an estimated measurement uncertainty of ±2 °C. The hot-edge temperature is approximately 82 °C and the far edge is approximately 51 °C.

Analysis

A linear model is

where a is the intercept in degrees Celsius and b is the temperature gradient in degrees Celsius per metre. A least-squares fit chooses a and b to minimize

This model is descriptive. It is not automatically a solution of the heat equation because boundary conditions, heat generation, contact resistance, and transient effects have not been established.

import numpy as np
import matplotlib.pyplot as plt

x = np.array([0.00, 0.05, 0.10, 0.15, 0.20])
T = np.array([82., 74., 66., 58., 51.])
u = np.full_like(T, 2.0)

coeff, cov = np.polyfit(x, T, 1, cov=True)
b, a = coeff
x_fit = np.linspace(x.min(), x.max(), 200)
T_fit = a + b*x_fit
residuals = T - (a + b*x)

The endpoint slope is (51 − 82)/(0.20 − 0) = −155 °C/m. Least-squares fitting gives a slope of approximately −154 °C/m and an intercept of approximately 81.8 °C; the fitted line is therefore close to the endpoint estimate. The residuals should be plotted and inspected rather than hidden. A systematic curved pattern would indicate that a linear model is inadequate.

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True,
                               figsize=(6.5, 6.0),
                               constrained_layout=True)
ax1.errorbar(x, T, yerr=u, fmt='o', capsize=3,
             color='black', label='Measurements')
ax1.plot(x_fit, T_fit, color='tab:red', label='Linear fit')
ax1.set_ylabel('Temperature (°C)')
ax1.legend(frameon=False)
ax1.grid(alpha=0.3)

ax2.axhline(0, color='black', linewidth=0.8)
ax2.scatter(x, residuals, color='tab:blue')
ax2.set_xlabel('Distance from heated edge (m)')
ax2.set_ylabel('Residual (°C)')
ax2.grid(alpha=0.3)
fig.savefig('plate_temperature.pdf')

Outcome and engineering interpretation

The data show a strong decrease in temperature with distance. The approximate gradient is −154 °C/m, but the uncertainty of ±2 °C means that the gradient should not be reported with many unjustified digits. The fitted line provides a useful first-order description over the tested interval. It does not establish the thermal conductivity of aluminium, because conductivity would require a heat flux, geometry, boundary conditions, and an appropriate heat-transfer model.

If the residuals are randomly distributed and mostly smaller than the measurement uncertainty, the linear approximation may be adequate for preliminary design. If residuals show a systematic trend, the team should investigate nonlinear conduction, heat loss to the surroundings, sensor placement, or insufficient equilibration time.

Learning points

The case illustrates five professional practices: units were explicit; uncertainty was displayed; a model was distinguished from raw data; residuals were examined; and engineering conclusions were limited to what the evidence supports.

7. Reproducibility, Ethics, and Professional Practice

A scientific script should be reproducible by another engineer. Record the Python version, package versions, input files, unit conventions, random seeds where relevant, and figure-generation commands. Use version control and meaningful filenames. Avoid manually editing plotted values after computation.

Data integrity is an ethical requirement. Do not delete inconvenient observations without a documented technical reason. Do not truncate axes to exaggerate a difference, use inconsistent scales between comparable panels, or omit uncertainty because it weakens a preferred conclusion. When preprocessing is performed—such as filtering, interpolation, or removal of outliers—state the method and its rationale.

AI-assisted programming can help generate boilerplate code, explain error messages, or suggest visualization styles. It can also produce incorrect formulas, silently change units, use deprecated functions, or invent citations. Engineers remain responsible for testing code, checking dimensions, validating numerical results, and protecting confidential data. Generated code should be treated as an unverified draft.

8. Limitations and Best Practices

Numerical output is not automatically accurate. Important limitations include floating-point round-off, truncation error, ill-conditioned matrices, insufficient sampling, interpolation bias, measurement noise, and model inadequacy. A visually smooth curve may be a consequence of interpolation rather than physical evidence.

Recommended practices are:

• Use SI units internally and label converted units clearly.
• Check array shapes before matrix operations.
• Prefer solvers such as np.linalg.solve to explicit matrix inversion.
• Test known limiting cases and simple benchmark problems.
• Compare numerical derivatives and integrals against analytical results when possible.
• Plot raw data as well as fitted or filtered data.
• Define error bars and uncertainty bands in captions.
• Use colour maps and scales appropriate to the data type.
• Preserve the code and data that produced every reported figure.
• Report only the precision justified by measurement and numerical uncertainty.

Chapter Summary

Scientific programming represents engineering quantities as scalars, vectors, matrices, and multidimensional arrays. NumPy provides efficient elementwise operations, matrix multiplication, linear-system solvers, norms, and numerical tools. Correct dimensions, units, conditioning, residuals, and model assumptions must be checked before interpreting results.

Matplotlib supports line, scatter, contour, heat-map, vector-field, and three-dimensional visualizations. A professional graph uses clear labels, units, suitable scales, accessible colours, meaningful legends, and an informative caption. Visualization is an analytical method: residual plots, uncertainty bands, and multidimensional maps can reveal behaviour that summary numbers conceal.

Uncertainty must be defined and communicated honestly. Statistical variation, measurement uncertainty, model inadequacy, and numerical error are distinct sources of limitation. Reproducible code, documented assumptions, ethical data handling, and critical evaluation are essential elements of engineering practice.

Key Terms

Array: An indexed collection of numerical values, often represented by a NumPy ndarray.

Broadcasting: NumPy's rule-based method for operating on arrays with compatible shapes.

Condition number: A measure of how sensitive a numerical problem is to perturbations in its inputs.

Contour plot: A graph showing curves of constant value over a two-dimensional domain.

Covariance: A measure of how two variables vary together.

Elementwise operation: An operation applied independently to corresponding array elements.

Ill-conditioned system: A system in which small input changes can cause large output changes.

Interpolation: Estimation of values between known data points.

Least squares: A fitting method that minimizes the sum of squared residuals.

Mesh grid: Coordinate arrays used to evaluate and visualize a function over a rectangular domain.

Matrix: A rectangular arrangement of numerical elements.

Norm: A mathematical measure of vector magnitude.

NumPy: A Python library for numerical arrays and scientific computation.

Residual: The difference between an observed or specified value and a model value.

Scalar: A single numerical quantity.

Standard error: An estimate of the uncertainty in a sample mean or another estimator.

Uncertainty interval: A stated range representing variability or uncertainty under specified assumptions.

Vector: An ordered collection of components with magnitude and, in many applications, direction.

Vectorization: Performing array operations without explicit element-by-element Python loops.

Visualization: Graphical representation of data or model results for analysis and communication.

Review Questions

Short Answer Questions

1. What is the difference between elementwise multiplication and matrix multiplication in NumPy?
2. Why is np.linalg.solve(A, b) generally preferred to explicitly computing A⁻¹b?
3. What information should be included in the label of an engineering graph?
4. Distinguish standard deviation from standard error of the mean.
5. Why should residuals be plotted after fitting a model?

Long Answer Questions

1. Explain how vectors and matrices support mechanical-engineering calculations, using one physical example for each.
2. Describe a complete workflow for converting experimental measurements into a reproducible scientific figure.
3. Compare line plots, scatter plots, contour plots, and heat maps, including an appropriate engineering use for each.
4. Explain the assumptions and limitations of first-order uncertainty propagation.
5. Discuss how an ill-conditioned matrix, inconsistent units, and insufficient sampling can each produce misleading numerical results.

Analytical / Application Questions

1. A displacement sensor records position at irregular time intervals. Design a Python-based procedure to estimate velocity, identify the effects of noise, and communicate uncertainty in the result.
2. A two-dimensional simulation produces pressure values on a rectangular grid. Propose a visualization strategy that reveals both absolute pressure and deviation from a reference pressure, and justify the colour maps and scales you would use.
3. An experimental temperature profile appears linear, but its residuals form a systematic curve. Explain at least three possible causes and propose additional analyses or measurements to distinguish among them.

Multiple Choice Questions with Answers

Question 1

Which NumPy expression performs matrix multiplication?

A. A * B

B. A @ B

C. A + B

D. A / B

Correct Answer: B. A @ B

Question 2

For the system A x = b, which condition is required for np.linalg.solve to produce a unique solution when A is square?

A. A must contain only integers

B. A must be diagonal

C. A must be nonsingular

D. b must be a matrix with two columns

Correct Answer: C. A must be nonsingular

Question 3

What does a residual represent in a model fit?

A. The model's independent variable

B. The difference between an observed value and a model value

C. The average of all measured values

D. The matrix condition number

Correct Answer: B. The difference between an observed value and a model value

Question 4

Which plot is most appropriate for showing a scalar field over two spatial coordinates?

A. Contour plot

B. Pie chart

C. Categorical bar chart

D. One-dimensional histogram only

Correct Answer: A. Contour plot

Question 5

In first-order uncertainty propagation, independent input variables have which covariance contribution?

A. It is always equal to one

B. It is omitted because covariance is zero

C. It is doubled

D. It replaces all derivative terms

Correct Answer: B. It is omitted because covariance is zero

Question 6

What is the primary purpose of a colour bar on a contour or pseudocolour plot?

A. To increase the number of data points

B. To identify the computer operating system

C. To map colours to quantitative values

D. To remove measurement uncertainty

Correct Answer: C. To map colours to quantitative values

Question 7

A Python list [1, 2, 3] multiplied by 2 produces which result?

A. [2, 4, 6]

B. [1, 2, 3, 2]

C. [1, 2, 3, 1, 2, 3]

D. [0.5, 1, 1.5]

Correct Answer: C. [1, 2, 3, 1, 2, 3]

Question 8

For a vector x, what does np.linalg.norm(x) calculate by default?

A. The Euclidean norm

B. The largest matrix eigenvalue only

C. The number of array dimensions

D. The arithmetic mean

Correct Answer: A. The Euclidean norm

Question 9

Which practice best supports reproducibility?

A. Manually editing plotted values after computation

B. Saving code, input data, package versions, and figure-generation settings

C. Removing all unexpected measurements without documentation

D. Using unlabeled axes to reduce visual clutter

Correct Answer: B. Saving code, input data, package versions, and figure-generation settings

Question 10

A temperature gradient calculated from 82 °C at 0 m and 51 °C at 0.20 m is approximately what value?

A. −6.2 °C/m

B. −31 °C/m

C. −155 °C/m

D. +155 °C/m

Correct Answer: C. −155 °C/m