Computational Methods
for Physicists and Materials Engineers
1
Basic programming I Python
two source codes file1.cpp and
file2.cpp in C++
generate object file file1.o
> g++ -c file1.cpp
generate object file file2.o
> g++ -c file2.cpp
link file1.o and file2.o to
generate executable prog.exe
> g++ -o prog.exe file1.o file2.o
run the executable file
> ./prog.exe
Programming process
Compiled languages (e.g. C, C++, Fortran) source code compiler object files
runtime library files
linker executable file run output
Programming process
Interpreted languages (e.g. python, JavaScript) source code
interpreter output
runtime library files
source code file.py in python
> python file.py
# head of file.py
# import library files
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack
Get started
# in “helloworld.py”
print(“Hello, World!”)
Indentation “tab”: a block of code
# define a function
Python syntax
def arithmetic_sum(term_start, comm_diff, num_terms): series_sum = 0
for i in range(num_terms):
series_sum += term_start + i * comm_diff
return series_sum
print(arithmetic_sum(4, 2, 3))
Comments
# this is a single-line comment
“””
this is a multi-line comment
another line
one more line
“””
Line continuation by “\”
def is_in_cube(point, xlo=-1, xhi=10, ylo=-20,\ yhi=3, zlo=0, zhi=100):
x, y, z = point
if xlo <= x <= xhi and ylo <= y <= yhi \ and zlo <= z <= zhi:
return "The point is located \
inside the cube."
else:
return "The point is outside the cube."
point = (0.1, 0.5, 30)
print(is_in_cube(point))
Multiple statements in one line
# separate statements by semicolon
a = 1; b = 2; print(a > b)
Python syntax
Import libraries
Python syntax
# at the head of each source code file
import numpy as np
import matplotlib.pyplot as plt
from domain import is_in_cube # import a function from another
# file “domain.py”
# in a separate file “domain.py”
def is_in_cube …
The libraries which we will use in this course: NumPy
Matrices https://numpy.org/ SciPy
Scientific computation https://www.scipy.org/ SymPy
Symbolic math https://www.sympy.org/en/index.html
matplotlib
plotting https://matplotlib.org/
a=str(3) #awillbe’3′ b=int(3) #bwillbe3 c = float(3) # c will be 3.0
Get the types
print(type(x))
print(type(y))
Variable name:
can only contain A‐z, 0‐9, and _
cannot start with a number
case‐sensitive (a and A are different variables)
Variables : Data types
Creating variables: no need to be declared with type x=3
y = “Hello”
Casting: specify data type of a variable
Category Numeric
Data type
Example
Boolean Text
bool
x = True; y = False
Sequence
list
tuple
range
x = [“a”, “b”, “c”]
x = (“a”, “b”, “c”)
x = range(6)
Set Mapping NumPy
set
dict
numpy.ndarray
x = {“apple”, “cherry”, “orange”}
x = {“element” : “Al”, “weight” : 27}
x = numpy.array([1, 2, 3])
int
float
complex
x = -365
x = 36.5; y = 5e2
x = 1 + 3j
str
x = “Hello World!”
Variables : Data types
Built‐in data types
representation of data + operations on data (methods)
Representation
x=3
y = -0.55 z=5
Methods (built‐in)
print(x + y)
print(x – y)
print(x * y)
print(x / y)
print(x // y)
print(x % y)
print(x ** y)
# sum
# difference
# product
# quotient; return float
# floor of quotient x / y
# remainder of quotient x / y
# x to the power y
Functions (built‐in)
print(abs(y))
print(int(y))
print(float(x))
# absolute value
# convert y to integer
# convert x to float
Variables : Data types and methods
Integers int & Floating points float representation of data + operations on data (methods)
Representation
x=3
y = -0.55 z=5
Functions (math Library) import math
print(math.sqrt(x)) print(math.log(x,z)) print(math.sin(x)) print(math.asin(y)) print(math.erf(y)) print(math.gamma(y))
# square root #basexlogarithm;bydefault,baseise # sine
# arc sine; note the allowed domain
# error function
# gamma function
Variables : Data types and methods
Integers int & Floating points float representation of data + operations on data (methods)
https://docs.python.org/3/library/math.html
Representation
x = “Hello World!” y=”ByeWorld ”
Methods (built‐in)
print(x[0] + x[6])
print(x[2:7])
# strings are arrays
# slicing
Array index
Variables : Data types and methods
Strings str
x Hello□World!
index 0 1 2 3 4 5 6 7 8 9 10 11 12 -12-11-10-9 -8 -7 -6 -5 -4 -3 -2 -1 0
Representation
x = “Hello World!” y=”ByeWorld ”
Methods (built‐in)
print(x[0] + x[6])
print(x[2:7])
# strings are arrays
# slicing
Array index
Variables : Data types and methods
Strings str
x Hello□World! index 0 1 2 3 4 5 6 7 8 9 10 11 12
-12-11-10-9 -8 -7 -6 -5 -4 -3 -2 -1 0
x[0] x[6] x[-12] x[-6]
Representation
x = “Hello World!” y=”ByeWorld ”
Methods (built‐in)
print(x[0] + x[6])
print(x[2:7])
# strings are arrays
# slicing
Array index
Variables : Data types and methods
Strings str
x Hello□World! index 0 1 2 3 4 5 6 7 8 9 10 11 12
-12-11-10-9 -8 -7 -6 -5 -4 -3 -2 -1 0
x[2:7]
x[-10:-5]
Representation
x = “Hello World!” y=”ByeWorld ”
Methods (built‐in)
print(x.upper())
print(x.lower())
print(y.strip())
print(x.split(” “))
print(x + ” then” + y)
# return string in upper case
# return string in lower case
# remove head and tail whitespaces
# split string by space
# string concatenation
Functions (built‐in)
print(len(x))
# number of items
String format
Variables : Data types and methods
Strings str
degree = 3
order = 4
value = 1.60217662
print(f”This is a degree-{degree} order-{order} ordinary \ differential equation.\n”) # \n means new line
print(f”The value is {value:0.3f}.”)
https://docs.python.org/3/library/string.html
Representation
Variables : Data types and methods
Array of data: list, tuple, dict, set Array of data (i): Lists
Ordered, changeable, and duplicates allowed
x = [“Gauss”, “Euler”, “Cauchy”, “Euler”] z = [3, 100, 24, 34, 1134]
y = [“abc”, 23.5, True, “Hello World!”, z]
Methods (built‐in)
print(y[3])
x[0] = “Einstein”
print(x[:2]) x.append(“Newton”); print(x) x.insert(3, “Dirac”); print(x) x.extend(x); print(x) x.remove(“Dirac”); print(x) print(x.pop(0)); print(x) z.reverse(); print(z) z.sort(reverse = True); print(z) a = z.copy(); print(a)
b = y + z; print(b)
# take one item
# replace item [0] by “Einstein”
# take a range of items
# append an item to the end
# insert Newton as the item [3]
# extend x with x
# remove the 1st “Dirac”
# remove item [0]
# reverse order
# sort descending
# copy z to a
# join x and y to form b
https://www.w3schools.com/python/python_lists_methods.asp
Nested lists: matrix
A = [[1, 2, 3],
[4, 5, 6],
[7, 8, 9]]
print(A[1][2])
# access an element of matrix
Variables : Data types and methods
Array of data: list, tuple, dict, set Array of data (i): Lists
Ordered, changeable, and duplicates allowed
[[A,A,A], 00 01 02
[A,A,A], 10 11 12
[A,A,A]] 20 21 22
Later we will show better representation of matrix: numpy.array
Representation
Methods (built‐in)
Variables : Data types and methods
Array of data (ii): Tuples
Ordered, unchangeable, and duplicates allowed
x = (“Gauss”, “Euler”, “Cauchy”, “Euler”) z = (3, 100, 24, 34, 1134)
y = (“abc”, 23.5, True, “Hello World!”, z)
print(y[3]) # take one item
print(x[:2]) # take a range of items
h, k, *L = z; print(h); print(k); print(L) # unpack a tuple
b = y + z; print(b)
print(x.count(“Euler”))
print(x.index(“Euler”))
# join y and z to form b
# number of occurrences of “Euler”
# get index of 1st “Euler”
cannot replace, (append, insert, extend), (remove, pop), (reverse, sort), copy
Representation
Methods (built‐in)
x.add(“Einstein”); print(x)
x.update(z); print(x)
x.remove(“Euler”); print(x)
print(x.pop()); print(x)
# similar to append
# similar to extend
# remove “Euler”
# remove a random item
# copy z to a
a = z.copy(); print(a)
b = y.union(z); print(b) c=y.intersection(z);print(c)
# union set of y and z #intersectionsetofyandz
Variables : Data types and methods
Array of data (iii): Sets
Unordered, unchangeable, and duplicates not allowed
x = {“Gauss”, “Euler”, “Cauchy”, “Euler”}
y = {“abc”, 23.5, True, “Hello World!”, (1, 2, 3)} # only tuple type of arrays can be an item of set z = {3, 100, 24, 34, 1134}
cannot access items by indices, replace, (reverse, sort), count, index https://www.w3schools.com/python/python_sets_methods.asp
Representation
x = { # key
“element”
: value
}
“structure”
Methods (built‐in)
print(x[“element”]) print(x.keys()) print(x.values()) print(x.items())
x[“weight”] = 27; print(x) print(x.pop(“weight”)); print(x) a = x.copy(); print(a)
# access value by key
# display all keys
# display all values
# display all items
# add item
Variables : Data types and methods
Array of data (iv): Dictionaries Unordered, changeable, and duplicates not allowed
: “fcc”,
“melting point” : 933.47,
: “Al”, Don’t forget comma
# remove the item with the key
# copy x to a
https://www.w3schools.com/python/python_dictionaries_methods.asp
Nested dictionaries
student_info = {
“student1” : {
Variables : Data types and methods
Array of data (iv): Dictionaries Unordered, changeable, and duplicates not allowed
},
“student2” : {
},
“student3” : {
}, print(student_info[“student2”][“name”])
}
“name” : “Hilbert”,
“major” : “math”,
“year” : 1862,
“name” : “Einstein”,
“major” : “physics”,
“year” : 1879,
“name” : “Russell”,
“major” : “philosophy”,
“year” : 1872,
import numpy as np
Variables : Data types and methods
arr = np.array([[1, 2, 3, 4 ], [5, 6, 7, 8],
[9, 10, 11, 12]])
print(arr); print(type(arr))
Access elements
We will learn NumPy in the next lecture
print(arr[0, 1])
print(arr[0:-1, 1:3])
NumPy array numpy.ndarray
https://numpy.org/
NumPy is a library supporting large, multi‐dimensional arrays and matrices, along with a large collection of high‐level functions to operate on these arrays and matrices
Representation
a = [“Gauss”, “Euler”] b=a
c = a.copy()
print(a == b)
# == : value equality
print(a == c)
print(a is b)
# is : are they exactly the same object?
# check memory address by id(x)
print(a is c)
Variables : Data types and methods
Comparison and Boolean bool Comparison returns Boolean
<, <=, >, >=, ==, !=, is, isnot
Boolean operations take in Boolean and return Boolean
x = True
y = False
print(x or y)
print(x and y)
print(not x)
# if x is False, then y, else x
# if x is False, then x, else y
# if x is False, then True, else False
or, and, not
Define constants: capital letters
create “constant.py”
PI = 3.1415926536
GRAVITY = 9.8
create “main.py”
import constant
print(constant.PI)
print(constant.GRAVITY)
Pre‐define constants: math library, scipy library import math
import scipy.constants
print(math.e)
print(math.pi)
print(math.inf)
print(scipy.constants.e)
print(scipy.constants.pi)
print(scipy.constants.k)
print(scipy.constants.hbar)
# Euler’s constant
# Pi
# +infinity
# Euler’s constant
# Pi
Variables : Constants
# Boltzmann constant
# Planck constant over 2*pi
https://docs.scipy.org/doc/scipy/reference/constants.html
Iterating over an array (string, list, tuple, set, dictionary, range)
Looping through string
x = “Hello World!”
for letter in x:
print(letter)
Looping through list, tuple, and set
list_arr = [“H”, “e”, “l”, “l”, “o”, ” “] tuple_arr = (“W”, “o”, “r”, “l”, “d”, “!”) for letter in list_arr:
print(letter, end=”)
for letter in tuple_arr:
print(letter, end=”)
print()
# loop over index
for index, _ in enumerate(list_arr): print(list_arr[index], end=”)
for index, _ in enumerate(tuple_arr): print(tuple_arr[index], end=”)
Looping
for
Looping through dictionary
material = { # key
“element”
: value
}
for key in material:
# loop over keys
# loop over values
# loop over items
“structure”
“melting point” : 933.47,
print(f”{key} : {material[key]}”) for value in material.values():
print(value)
for key, value in material.items():
print(f”{key} : {value}”)
Looping through range
range(start, stop, step)
An array of integers (by default start = 0):
start, start + step, start + 2*step, ···, step – 1
for x in range(3, 20, 2):
print(x)
Looping
for
: “Al”,
: “fcc”,
2 21 1
22 2
2N N
Looping
for
Example: multiplication of a matrix with a vector
A AAxb 11 12 1N 1 1
21 22 2N 2 2 A AAxb
A A Axb
M1 M2 MN N M b A x A x A x
1 11 1 12 2 1N N b A x A x A x
b AxAxAx j1
i
i1 1 i2 2 iN N
b M
A xA xA x M1 1 M2 2 MN N
N b Ax
i
(i 1,2,,M)
ij j
import sys
A = [[1, 2], [3, 4], [5, 6]]
x = [7, 8]
num_row = len(A)
num_col = len(A[0])
if len(x) != num_col:
sys.exit(“Error: dimensions do not match!”)
Looping
b = []
for i in range(num_row):
b_sum = 0
for j in range(num_col):
N
b Ax
b_sum = A[i][j] * x[j]
b.append(b_sum)
ij j
(i 1,2,,M)
print(b)
for
Example: multiplication of a matrix with a vector
i
j1
import sys
import numpy as np
A = np.array([[1, 2], [3, 4], [5, 6]])
x = np.array([[7], [8]])
#———————————————–
num_row, num_col = A.shape
if x.shape[0] != num_col:
sys.exit(“Error: dimensions do not match!”) b = np.zeros(shape=(num_row, 1))
Looping
for i in range(num_row):
for j in range(num_col):
b[i] += A[i, j] * x[j]
N
b Ax
#———————————————– # the above part is equivalent to
# b = np.matmul(A, x)
ij j
(i 1,2,,M)
print(b)
for
Example: multiplication of a matrix with a vector
i
j1
Execute statements as long as a condition is true
# calculate factorial N!
N=6
i=1 factorial = 1 while i <= N:
factorial *= i
i += 1
print(factorial)
trajectory = ['left', 'up', 'down', 'left', 'right'] while trajectory:
print(trajectory.pop(0))
# a list is False when it is empty
Looping
while
import numpy as np
def step_func(x):
# define a Heaviside step function
if x < 0:
return 0
elif x > 0:
return 1
else:
return 0.5
dx =0.2
x = np.arange(-1, 1 + dx, dx) y = x.copy()
for i in range(len(x)):
y[i] = step_func(x[i])
print(y)
Condition
if, elif, else
import numpy as np
# have defined step_func(x)
def step_func_arr(x):
# allow element-wise operation
# x can be a number, a list or a numpy.ndarray if type(x) is list or type(x) is np.ndarray:
y = x.copy()
for i in range(len(x)):
y[i] = step_func(x[i])
return y
else:
return step_func(x)
dx = 0.2
x = np.arange(-1, 1 + dx, dx)
y = step_func_arr(x)
print(y)
Condition
if, elif, else
Stop a loop
def sawtooth(x):
sum = 0
xmin = 0.01; xmax = 0.99
fmin = sawtooth(xmin); fmax = sawtooth(xmax) tolerance = 1e-15; step = 0; print_step = 5 whileTrue: #mustcontainbreak
x = (xmin + xmax) / 2
f = sawtooth(x)
if f * fmin > 0:
xmin = x; fmin = sawtooth(xmin)
else:
xmax = x; fmax = sawtooth(xmax)
step += 1
Branching
if step % print_step == 0:
print(f”Step {step}: priori error {x – 0.5}.”)
if xmax – xmin < tolerance:
print(f"The approximate zero is {x}."); break
break
for n in range(1, 1000000):
sum -= math.sin(2 * math.pi * n * x) / (math.pi * n)
return sum
Skip a cycle
# print all prime numbers in a range
# prime number must be greater than 1
if num <= 1:
continue
for i in range(2, num):
if num % i == 0:
else:
print(num)
break
Branching
lower = -10
upper = 100
print(f"Prime numbers between {lower} and {upper} are: ")
for num in range(lower, upper + 1):
continue
writing this operation as a function
for num in range(lower, upper + 1):
# prime number must be greater than 1
if num <= 1:
continue
for i in range(2, num):
if num % i == 0:
Functions
When we want to repeatedly do an operation, we should consider
# define a function to print all prime numbers in a range
def prime_print(lower, upper):
break
else:# this "else" goes with "for", not "if"
print(num)
lower = -10; upper = 100
print(f"Prime numbers between {lower} and {upper} are: ") prime_print(lower, upper)
lower = 200; upper = 300
print(f"Prime numbers between {lower} and {upper} are: ") prime_print(lower, upper)
def
A function with both input (arguments) and output (returns)
# define a function to tell if a number is prime
def is_prime(num):
# prime number must be greater than 1
if num <= 1:
return f"{num} is not a prime."
for i in range(2, num):
if num % i == 0:
return f"{num} is not a prime."
return f"{num} is a prime."
nums = [-12, 0, 11, 311, 2445]
for num in nums:
print(is_prime(num))
Functions
def
# in a file "number_tool.py"
def is_prime(num):
# prime number must be greater than 1
if num <= 1:
return f"{num} is not a prime."
for i in range(2, num):
if num % i == 0:
return f"{num} is not a prime."
return f"{num} is a prime."
Using the module
# in the main code "main.py"
import number_tool
nums = [-12, 0, 11, 311, 2445]
for num in nums:
print(number_tool.is_prime(num))
Modules
A separate file containing a set of variables or functions you can import in the main code
Creating a module
Creating a module
# in a file "constants.py"
material_constants = {
"Al" : {
"crystal_structure" : "Fcc", "lattice_constant" : 4.05, },
}
Using the module
# in the main code "main.py"
from constants import material_constants as matconst
for elem, property in matconst.items():
structure = property["crystal_structure"].lower() lattconst = property["lattice_constant"]
if structure == "bcc":
atom_vol = lattconst**3 / 2
elif structure == "fcc":
atom_vol = lattconst**3 / 4
else:
Modules
"W" : { "crystal_structure" : "bCc", "lattice_constant" : 3.16, },
print(f"{structure} is an unknown structure."); break print(f"The atomic volume of {elem} is {atom_vol} A^3.")
stocks = [] ‘read’
file = open("HSI.csv", 'r')
# read the file line by line
for line_index, line_str in enumerate(file):
if line_index == 0:
continue # skip the 1st line
stocks.append(line_str)
file.close()
print(stocks)
File handling
Open and read a file
In the current folder, there is a file which contains data: HSI.csv
https://finance.yahoo.com/ search Hang Seng index
Each line is read in form of a string
E.g. the 1st item of stocks is "2001-01-01,15089.849609,16255.110352,
14512.709961,16102.349609,16102.349609,0\n"
But the data we want are
"2001-01-01" (str), 15089.849609 (float), 16255.110352 (float), 14512.709961 (float), 16102.349609 (float), 16102.349609 (float)
How to split the string to a set of data?
Open and read a file How to split the string to a set of data?
string = string.replace(' ', ' ')
File handling
def split_string_to_data(string):
string = string.replace('\n', '') # delete tail '\n' string = string.replace(',', ' ') # replace ',' by ' ' string = string.replace('\t', ' ') # replace '\t' by ' ' string = string.replace(';', ' ') # replace ';' by ' ' while ' ' in string:
# replace multiple spaces by one space
# split with delimiter ' ' and store them in a list
var_list = string.split(' ')
while '' in var_list:
# remove empty strings from the list
var_list.remove('')
return var_list
stocks = {
"index" : [],
"date" : [],
"open" : [],
"high" : [],
"low" : [],
"close" : [],
} # the dictionary is a container for data
file = open("HSI.csv", 'r')
for line_index, line_str in enumerate(file): if line_index == 0:
continue # skip the 1st line
line_list = split_string_to_data(line_str) stocks["index"].append(line_index - 1) stocks["date"].append(line_list[0]) stocks["open"].append(float(line_list[1])) stocks["high"].append(float(line_list[2])) stocks["low"].append(float(line_list[3])) stocks["close"].append(float(line_list[4]))
file.close()
File handling
Open and read a file
If the file exists, overwrite; if not, create one
file = open("HSI_copy.csv", 'w')
file.write(f"Date, Open, High, Low, Close\n") for i in stocks["index"]:
File handling
Open and write a file
date = stocks["date"][i]
open = stocks["open"][i]
high = stocks["high"][i]
low = stocks["low"][i]
close = stocks["close"][i]
file.write(f"{date}, {open:>12.6f}, {high:>12.6f}, {low:>12.6f}, {close:>12.6f}\n”)
file.close()
{ open : > 12 . 6 f }
right justify 12 digits fixed point
6 digits after decimal point
The Zen of Python 襌 import this
Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren’t special enough to break the rules. Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one – and preferably only one – obvious way to do it. Although that way may not be obvious at first unless you’re Dutch. Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it’s a bad idea.
If the implementation is easy to explain, it may be a good idea. Namespaces are one honking great idea – let’s do more of those!