MATH 446/546: Scientific Computing in Python Lecture 0

MATH 446/546: Scientific Computing in Python

Python Fundamentals for Scientific Computing
Lecture 0, Fall 2026

Overview
This lecture note introduces the basic Python programming skills needed for this course. We assume
no prior experience with Python. Our goal is not to provide a comprehensive introduction to Python
as a general-purpose language, but rather to cover the essential ingredients for implementing and
exploring the numerical methods we will study throughout the semester. In particular, we focus
on the NumPy library, which will be our primary tool for array-based computation throughout the
course.
A companion Jupyter notebook (Lecture [Link]) accompanies this note. All code examples
below appear in the notebook, where you can run, modify, and experiment with them. You are
encouraged to work through the notebook interactively after reading this material.

1 Variables and Arithmetic

In Python, we assign values to variables using the = sign. Unlike many compiled languages, Python
does not require us to declare variable types in advance.
1 x = 3.0
2 y = 2
3 z = x + y
4 print ( z ) # Output : 5.0

1.1 Arithmetic operators

Python supports the standard arithmetic operations:

Operator Description Example

+ Addition 3 + 2→5
- Subtraction 3 - 2→1
* Multiplication 3 * 2→6
/ Division 7 / 2 → 3.5
** Exponentiation 3**2 → 9
% Modulo (remainder) 7 % 2→1
// Integer division 7 // 2 → 3

1
MATH 446/546: Scientific Computing in Python Lecture 0

1.2 Numeric types: int and float

Python distinguishes between integers (int) and floating-point numbers (float). This distinction
will become important when we study floating-point arithmetic later in the course.
1 a = 4 # int
2 b = 4.0 # float
3 print ( type ( a ) ) # < class ’ int ’>
4 print ( type ( b ) ) # < class ’ float ’>

Note that dividing two integers with / always produces a float:

1 print (4 / 2) # Output : 2.0 ( float , not int )
2 print ( type (4 / 2) ) # < class ’ float ’>

Example 1.1. Consider the polynomial p(x) = x3 − 2x + 1. We can evaluate it at x = 1.5 as

follows:
1 x = 1.5
2 p = x **3 - 2* x + 1
3 print ( p ) # Output : 1.375

We can verify this by hand: 1.53 − 2(1.5) + 1 = 3.375 − 3 + 1 = 1.375.

1.3 Boolean values and comparisons

Python has a Boolean type bool with values True and False. Comparison operators produce
Boolean values:
1 print (3 > 2) # True
2 print (3 == 2) # False
3 print (3 != 2) # True
4 print (3 >= 3) # True

Remark. A single = is assignment, while a double == is comparison. Writing x = 3 sets x to 3;

writing x == 3 tests whether x equals 3. This is a common source of bugs for beginners.

2 Loops and Conditionals

2.1 The for loop
A for loop iterates over a sequence. The most common use in scientific computing is to repeat a
computation a fixed number of times. In Python, this is done with the range() function.
1 for i in range (5) :
2 print ( i )
3 # Output : 0 , 1 , 2 , 3 , 4

Note that range(n) produces integers from 0 to n − 1. More generally, range(a, b) produces
integers from a to b − 1.
Remark. Python uses indentation (typically 4 spaces) to define blocks of code. All lines inside a
loop or conditional must be indented to the same level. This replaces the curly braces used in
languages like C or Java.

2
MATH 446/546: Scientific Computing in Python Lecture 0

PN 1
Example 2.1. We can compute the partial sum SN = k=1 k2 as follows:
1 N = 1000
2 S = 0.0
3 for k in range (1 , N + 1) :
4 S = S + 1.0 / k **2
5 print ( S ) # Output : 1.6 43 93 45 66 68 15 61 5

The exact value of the infinite series is π 2 /6 ≈ 1.6449340668 . . ., so with N = 1000 terms we have
roughly 3 correct digits.

2.2 The while loop

A while loop repeats as long as a condition is true. This is particularly natural for iterative
algorithms, where we want to repeat until some convergence criterion is met.
√
Example 2.2 (The Babylonian method for a). One of the oldest known algorithms computes
√
a by the iteration  
1 a
xn+1 = xn + , n = 0, 1, 2, . . .
2 xn
starting from an initial guess x0 > 0. This is known as the Babylonian method (or Heron’s method).
We will later see that this is a special case of Newton’s method. For now, we implement it using a
while loop:
1 a = 2.0
2 x = 1.0 # initial guess
3 tol = 1e -10 # tolerance
4
5 while True :
6 x_new = 0.5 * ( x + a / x )
7 if abs ( x_new - x ) < tol :
8 break
9 x = x_new
10
11 print ( x_new ) # Output : 1.4 14 21 35 62 37 30 95 1
√
The result agrees with 2 = 1.41421356237 . . . to machine precision.

2.3 Conditionals: if, elif, else

Conditional statements allow us to execute different code depending on whether a condition is true
or false.
1 x = -3.5
2
3 if x > 0:
4 print ( " x is positive " )
5 elif x == 0:
6 print ( " x is zero " )
7 else :
8 print ( " x is negative " )
9 # Output : x is negative

3
MATH 446/546: Scientific Computing in Python Lecture 0

Example 2.3. The absolute value function can be implemented as:

1 def my_abs ( x ) :
2 if x >= 0:
3 return x
4 else :
5 return -x

(We will discuss functions in detail in the next section.)

3 Functions
In this course, we will frequently package algorithms as functions. A Python function is defined
with the def keyword:
1 def function_name ( argument1 , argument2 ) :
2 # body of the function
3 result = ...
4 return result

The return statement specifies the output of the function.

Example 3.1 (Evaluating a polynomial). We can turn the polynomial evaluation from Example 1.1
into a reusable function:
1 def p ( x ) :
2 return x **3 - 2* x + 1
3
4 print ( p (0) ) # Output : 1
5 print ( p (1) ) # Output : 0
6 print ( p (1.5) ) # Output : 1.375

Example 3.2 (The Babylonian method as a function). We can package the Babylonian method
from Example 2.2 into a function that takes the input a and a tolerance:
1 def my_sqrt (a , tol =1 e -10) :
2 """
3 Compute sqrt ( a ) using the Babylonian method .
4
5 Parameters :
6 a : a positive number
7 tol : stopping tolerance ( default : 1e -10)
8
9 Returns :
10 Approximation of sqrt ( a )
11 """
12 x = a # initial guess
13 while True :
14 x_new = 0.5 * ( x + a / x )
15 if abs ( x_new - x ) < tol :
16 return x_new
17 x = x_new
18
19 print ( my_sqrt (2) ) # 1.4 14 21 35 62 37 30 95 1

4
MATH 446/546: Scientific Computing in Python Lecture 0

20 print ( my_sqrt (9) ) # 3.0

21 print ( my_sqrt (2 , tol =1 e -4) ) # less accurate

Several things to note about functions:

• The triple-quoted string """...""" at the beginning of the function body is called a docstring.
It documents what the function does. Writing docstrings is good practice and will be expected
in this course.
• The parameter tol=1e-10 has a default value. If the caller does not supply a value for tol,
Python uses 1e-10. This is convenient for parameters that usually take a standard value.
• A function can take any number of arguments and can return any type of object.

4 NumPy
Python’s built-in arithmetic operates on single numbers. For scientific computing, we almost always
work with arrays of numbers—vectors, matrices, grids of function values—and we need to perform
operations on entire arrays efficiently. The NumPy library provides exactly this capability. It will
be our most important computational tool throughout the course.
We import NumPy using the standard abbreviation:
1 import numpy as np

4.1 Creating arrays

A NumPy array is an ordered collection of numbers of the same type, analogous to a vector in
mathematics. There are several common ways to create arrays.
From a Python list. The most direct method is to pass a list of numbers to [Link]:
1 v = np . array ([1.0 , 4.0 , 9.0 , 16.0])
2 print ( v ) # [ 1. 4. 9. 16.]
3 print ( len ( v ) ) # 4

Evenly spaced points. The function [Link](a, b, n) creates an array of n evenly spaced
points from a to b, inclusive:
1 x = np . linspace (0 , 1 , 5)
2 print ( x ) # [0. 0.25 0.5 0.75 1. ]

This is one of the most frequently used functions in this course: whenever we need to evaluate a
function on an interval, we typically begin by creating a grid with [Link].
The function [Link]. Similar to Python’s range, the function [Link](a, b, h) creates
an array from a to b (exclusive) with step size h:
1 x = np . arange (0 , 1 , 0.2)
2 print ( x ) # [0. 0.2 0.4 0.6 0.8]

Special arrays. NumPy provides functions for creating arrays of zeros, ones, or uninitialized
entries:

5
MATH 446/546: Scientific Computing in Python Lecture 0

1 z = np . zeros (5) # [0. 0. 0. 0. 0.]

2 w = np . ones (5) # [1. 1. 1. 1. 1.]
3 e = np . empty (5) # 5 - element array ( uninitialized )

The function [Link] is especially useful for pre-allocating an array that will be filled in by a
loop, as we will frequently do when implementing iterative methods.

4.2 Vectorized operations

The key advantage of NumPy arrays is vectorized operations: arithmetic is applied element-wise,
without the need for explicit loops.
1 x = np . array ([1.0 , 2.0 , 3.0 , 4.0])
2 print ( x + 10) # [11. 12. 13. 14.]
3 print ( x * 2) # [ 2. 4. 6. 8.]
4 print ( x **2) # [ 1. 4. 9. 16.]
5 print (1.0 / x ) # [1. 0.5 0.333 0.25 ]

Operations between two arrays of the same length are also element-wise:
1 x = np . array ([1.0 , 2.0 , 3.0])
2 y = np . array ([4.0 , 5.0 , 6.0])
3 print ( x + y ) # [5. 7. 9.]
4 print ( x * y ) # [ 4. 10. 18.]

Remark. Vectorized operations are not merely a notational convenience. NumPy executes these
operations in optimized compiled code (written in C), making them dramatically faster than equiv-
alent Python for loops. For large arrays, the difference can be a factor of 100 or more. Throughout
this course, we will prefer vectorized operations over explicit loops whenever possible.

4.3 Mathematical functions

NumPy provides standard mathematical functions that operate element-wise on arrays:
1 x = np . array ([0.0 , np . pi /6 , np . pi /4 , np . pi /3 , np . pi /2])
2 print ( np . sin ( x ) ) # [0. 0.5 0.707 0.866 1. ]
3 print ( np . cos ( x ) ) # [1. 0.866 0.707 0.5 0. ]
4 print ( np . exp ( x ) ) # [1. 1.69 2.19 2.85 4.81 ]

The most commonly used functions include:

NumPy function Mathematical equivalent

[Link](x), [Link](x), [Link](x) sin x, cos x, tan x
[Link](x), [Link](x) ex , ln x
√
[Link](x) x
[Link](x) |x|
[Link] π (constant)
np.e e (constant)

Remark. Python’s standard library includes a math module with functions like [Link], [Link],
etc. These only accept single numbers, not arrays. In this course, we will always use the NumPy
versions ([Link], [Link], etc.), which work on both single numbers and arrays.

6
MATH 446/546: Scientific Computing in Python Lecture 0

4.4 Indexing and slicing

Individual elements of an array are accessed by their index, starting from 0:
1 x = np . array ([10 , 20 , 30 , 40 , 50])
2 print ( x [0]) # 10 ( first element )
3 print ( x [2]) # 30 ( third element )
4 print ( x [ -1]) # 50 ( last element )

A slice extracts a subarray using the syntax x[start:stop], where the element at index stop
is not included :
1 print ( x [1:4]) # [20 30 40] ( indices 1 , 2 , 3)
2 print ( x [:3]) # [10 20 30] ( first 3 elements )
3 print ( x [2:]) # [30 40 50] ( from index 2 onward )

We can also modify individual entries or slices:

1 x [0] = 99
2 print ( x ) # [99 20 30 40 50]

Indexing and slicing will be essential when we implement numerical methods. For instance, an
iterative method that stores all its approximations x0 , x1 , . . . , xN in an array will use indexing to
access and update individual iterates.

4.5 Useful array operations

NumPy provides many functions for working with arrays. Here are some that we will use frequently:
1 x = np . array ([3.0 , 1.0 , 4.0 , 1.0 , 5.0])
2 print ( np . max ( x ) ) # 5.0
3 print ( np . min ( x ) ) # 1.0
4 print ( np . sum ( x ) ) # 14.0
5 print ( np . mean ( x ) ) # 2.8
6 print ( np . argmax ( x ) ) # 4 ( index of the maximum )

The function [Link]([Link](x)) computes the ℓ∞ norm ∥x∥∞ = maxi |xi |, which we will use
frequently to measure errors:
1 x = np . array ([0.1 , -0.5 , 0.3])
2 print ( np . max ( np . abs ( x ) ) ) # 0.5

Example 4.1. To evaluate p(x) = x3 − 2x + 1 at 100 evenly spaced points on [−2, 2]:
1 x = np . linspace ( -2 , 2 , 100)
2 y = x **3 - 2* x + 1

No loop is needed—NumPy handles the element-wise operations automatically. This is both cleaner
to write and significantly faster for large arrays than the equivalent for loop:
1 # Equivalent but slower ( avoid this style ) :
2 y = np . zeros (100)
3 for i in range (100) :
4 y [ i ] = x [ i ]**3 - 2* x [ i ] + 1

7
MATH 446/546: Scientific Computing in Python Lecture 0

4.6 Two-dimensional arrays (matrices)

NumPy also supports two-dimensional arrays, which represent matrices. We will need these when
we study systems of equations and interpolation.
1 A = np . array ([[1 , 2 , 3] ,
2 [4 , 5 , 6]])
3 print ( A . shape ) # (2 , 3) -- 2 rows , 3 columns
4 print ( A [0 , 1]) # 2 ( row 0 , column 1)
5 print ( A [1 , :]) # [4 5 6] ( entire row 1)
6 print ( A [: , 0]) # [1 4] ( entire column 0)

Special matrix constructors include:

1 Z = np . zeros ((3 , 3) ) # 3 x3 zero matrix
2 I = np . eye (3) # 3 x3 identity matrix
3 D = np . diag ([1 , 2 , 3]) # 3 x3 diagonal matrix

Matrix–vector multiplication is performed with the @ operator:

1 A = np . array ([[1 , 2] , [3 , 4]])
2 x = np . array ([1 , 1])
3 print ( A @ x ) # [3 , 7]

Remark. The * operator on two-dimensional arrays performs element-wise multiplication, not ma-
trix multiplication. Always use @ (or equivalently [Link]) for matrix–vector and matrix–matrix
products.

5 Matplotlib
Visualization is an essential part of scientific computing. The Matplotlib library provides plotting
capabilities. We import it as:
1 import matplotlib . pyplot as plt

The most common plot type is a line plot of y = f (x), created with [Link].

Example 5.1 (Plotting a polynomial and identifying roots). We plot p(x) = x3 − 2x + 1 on the
interval [−2, 2]:
1 x = np . linspace ( -2 , 2 , 200)
2 y = x **3 - 2* x + 1
3
4 plt . figure ( figsize =(6 , 4) )
5 plt . plot (x , y , ’b - ’ , linewidth =1.5 , label = r ’ $p ( x ) = x ^3 - 2 x + 1 $ ’)
6 plt . axhline ( y =0 , color = ’k ’ , linewidth =0.5)
7 plt . xlabel ( ’x ’)
8 plt . ylabel ( ’p ( x ) ’)
9 plt . title ( ’A cubic polynomial ’)
10 plt . legend ()
11 plt . grid ( True , alpha =0.3)
12 plt . tight_layout ()
13 plt . show ()

8
MATH 446/546: Scientific Computing in Python Lecture 0

The string ’b-’ specifies a blue solid line; other common choices include ’r--’ (red dashed), ’g:’
(green dotted), and ’ko’ (black circles). The prefix r before a string (as in r’$...$’) indicates a
raw string, which allows us to include LATEX formatting in labels and titles.
From the plot, we can visually identify that p(x) has roots near x ≈ −1.6, x = 1, and x ≈ 0.6.
In the next lecture, we will develop systematic numerical methods to find these roots precisely.
PN 1
Example 5.2 (Visualizing convergence). We return to the partial sums SN = k=1 k2 from
2
Example 2.1 and plot the error |SN − π /6| as a function of N :
1 N_values = np . arange (1 , 101)
2 S = 0.0
3 errors = np . zeros (100)
4 exact = np . pi **2 / 6
5
6 for i , N in enumerate ( N_values ) :
7 S = S + 1.0 / N **2
8 errors [ i ] = abs ( S - exact )
9
10 plt . figure ( figsize =(6 , 4) )
11 plt . semilogy ( N_values , errors , ’b - ’ , linewidth =1.5)
12 plt . xlabel ( ’N ’)
13 plt . ylabel ( ’ Error ’)
14 plt . title ( r ’ Convergence of $ \ sum_ { k =1}^{ N } 1/ k ^2 $ to $ \ pi ^2/6 $ ’)
15 plt . grid ( True , alpha =0.3)
16 plt . tight_layout ()
17 plt . show ()

The [Link] function uses a logarithmic scale on the y-axis, which is the standard way to
visualize convergence behavior in numerical analysis. On this scale, the error decreases roughly as
a straight line, suggesting that |SN − π 2 /6| ≈ C/N for some constant C. We will make this type
of analysis precise throughout the course.

Exercises
1. Evaluate the expression
310 − 210
310 + 210
in Python and verify your result by hand (or with a calculator).

2. Write a for loop that computes 10! (10 factorial). Verify your answer using the formula
10! = 3628800.

3. Modify the Babylonian method function my sqrt from Example 3.2 to additionally print the
√
iterate xn and the error |xn − a| at each step. (You may use [Link](a) to compute the
exact value for comparison.) Run it for a = 2 and observe how quickly the iterates converge.

4. Write a function horner(coeffs, x) that evaluates a polynomial using Horner’s method.

Given coefficients a0 , a1 , . . . , an (so that p(x) = an xn + · · · + a1 x + a0 ), Horner’s method
computes:


p(x) = a0 + x a1 + x a2 + · · · + x(an−1 + x · an ) · · · .
Test your function on p(x) = x3 − 2x + 1 at x = 1.5.

9
MATH 446/546: Scientific Computing in Python Lecture 0

5. Use NumPy and Matplotlib to plot the functions sin(x), cos(x), and sin(x) cos(x) on the
interval [0, 2π] in a single figure with a legend. Experiment with different line colors and
styles.

6. Plot the function

sin(x)
f (x) =
x
on the interval [−4π, 4π]. What happens at x = 0? How does Python handle this?
   
2 1 1
7. Let A = and b = . Use NumPy to:
1 3 2

(a) Create the matrix A and vector b as NumPy arrays.

(b) Compute Ab using the @ operator.
(c) Solve the linear system Ax = b using [Link](A, b) and verify that Ax = b
by computing the residual ∥Ax − b∥∞ .
1
8. Create a 10 × 10 matrix A with entries Aij = i+j−1 for i, j = 1, . . . , 10. (This is the Hilbert
matrix, which arises in polynomial fitting and is notoriously ill-conditioned.) Print the matrix
and compute its largest and smallest entries using [Link] and [Link].

What’s Next
In the next lecture, we will begin our study of solving nonlinear equations. Given a function f , we
seek values x∗ such that f (x∗ ) = 0. The plot in Example 5.1 already hints at the problem: we can
see approximately where the roots are, but how do we compute them accurately and efficiently?
We will study the bisection method and fixed-point iteration as our first systematic approaches.

10