by: Tomasz Romańczukiewicz
web: th.if.uj.edu.pl/~trom/
rm B-2-03
Programing languages:
It is usually listed as the second best tool for a given problem:
Task: Calculate a sume $s=\scriptsize\displaystyle\sum_{n=1}^Nn$ for some large $N$
%%time
N = int(1e7)
s = 0
for n in range(1,N+1):
s += n
import numpy as np
%%time
N = int(1e7)
s = np.sum(np.arange(1,N+1))
open iterative console: type python3 or ipython3 in a terminal)
1+2
0.1+0.2
0.1+0.2-0.3
type(1+2)
type(0.1+0.2)
1/2
type(1/2)
1//2
type(1//2)
2**5
6%4
help(print)
Iterative console is fun (automatic output) but for serious job we need a proper editor (IDE)
# Ordinary comment
"""
Multiline documentation
"""
print("Hello world") # print out a string
a = 5 # assign a variable
b = 6
print(a, b) # print out numbers
a, b = 3, 2 # multiple assignment, a tupple
print("a = ", a, " b = " , b)
a, b = b, a # simple swap
print("a = ", a, " b/a = ", b/a, "b//a = ", b//a) # unlike in C,Java etc b/a is not integer,
# but b//a is an integer
print("Formatted output: {:10d}, {:.4f}".format(a, b/a))
print("Conversion from floats to ints:", int(1000*a/b))
print("modulo division:", 15%6)
""" Logic operators """
to_be = True
print("1.", to_be or not to_be)
to_be = False
print("2.", to_be or not to_be)
print("3.", to_be and not to_be)
to_be = 3
print("4.", not to_be)
print("5.", 7 | 2, 7&2, ~7)
print("6.", 1<<3, 3<<3, 25>>2)
print("7.", (7%3!=0), (7%3!=0)and(7%2==0), (666<=777))
""" From 3.6 version """
# import math as mth # including the full math module with prefix mth
# print(f"sin(1) = {mth.sin(1)}")
# print(f"cos(1) = {mth.cos(1)}")
# import numpy as np # including the full numpy module with prefix np, in order not to confuse with math
# print(f"sin(1) = {np.sin(1)}")
# print(f"cos(1) = {np.cos(1)}")
# from math import exp, cosh # including only two functions but without prefix
# print(f"exp(1) = {exp(1)}")
# print(f"cosh(1) = {cosh(1)}")
import math as mth # including the full math module with prefix mth
print("sin(1) = ", mth.sin(1))
print("cos(1) = ", mth.cos(1))
import numpy as np # including the full numpy module with prefix np, in order not to confuse with math
print("sin(1) = ", np.sin(1))
print("cos(1) = ", np.cos(1))
from math import exp, cosh # including only two functions but without prefix
print("exp(1) = ", exp(1))
print("cosh(1) =", cosh(1))
Calculate the kinetic energy using relativistic and newtonian definitions
$$\scriptsize E_{rel} = \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^2$$$$\scriptsize E_{Newt} = \frac{1}{2}mv^2$$for $v=1 \frac{\text{m}}{\text{s}}$, $c=3\cdot 10^9 \frac{\text{m}}{\text{s}}$, $m=1$kg
from math import sqrt
v, c, m = 1, 3e9, 1
Erel = m*c**2/sqrt(1-v**2/c**2)-m*c**2
ENewt = 0.5*m*v**2
print("E_rel = ", Erel, "\nE_Newt = ", ENewt)
What's wrong?
v**2/c**2
1-v**2/c**2
Standard floats (64 bits) can handle only up to 16 digits So the accurate equation reduces to
$$\scriptsize E_{rel} = \frac{mc^2}{\sqrt{1-\frac{v^2}{c^2}}}-mc^2\to\frac{mc^2}{1}-mc^2=0$$which is painfully wrong.
Due to the cancellation errors sometimes it is better to use approximation rather than the full formula
But this approximation is wrong for large values of $v$.
Is there any way to have one formula which works in both cases?
With a bit of magic (short multiplication formulas)
$$\scriptsize E_{rel} = mc^2\left(\frac{c}{\sqrt{c^2-v^2}}-1\right) = mc^2 \frac{\frac{c^2}{{c^2-v^2}}-1}{\frac{c}{\sqrt{c^2-v^2}}+1} = \frac{mv^2}{\sqrt{1-\frac{v^2}{c^2}} + 1-\frac{v^2}{c^2}} $$g = 1-v**2/c**2
E3 = m*v**2/(sqrt(g)+g)
print(E3)
The above formula is mathematically equivalent with the definition but cancellation does not produce errors for small values of $v$.
""" Simple calculations: but what can go wrong? """
from numpy import cosh
print("1. ", 0.3+0.2-0.5) # this is 0
print("2. ", 0.1+0.2-0.3) # this is not 0, beware of round of errors
x = 20
print("3. ", (1-cosh(x))/(1+cosh(x)))
x = 800
print("4. ", (1-cosh(x))/(1+cosh(x))) # this produces overflow error
print("5. ", (1/cosh(x)-1)/(1/cosh(x)+1)) # this produces overflow error
""" List examples """
v = [3, 4, 5, 8, -1] # example array
print ("1:", v)
print ("2:", v[0], v[4], v[-2]) # elements: [(3), (4), 5, (8), -1]
# indexes: 0 1 2 3 4
# negative ind: -5 -4 -3 -2 -1
v.append(7) # adding a new element
print("3:", v)
v.insert(2, 6) # inserting at certain position
print("4:", v)
print("5:", v[2:4]) # slices [3, 4, (6, 5), 8, -1, 7]
v.sort() # hmm, what could that be?
print("6:", v)
w = list(range(1, 15, 3)) # each 3rd integer number from [1,15)
print (" 7:", w)
print(" 8:", w[2:5]+v[0:4:2]) # join list slices [1, 4, (7, 10, 13)] [(-1), 3, (4), 5, 6, 7, 8]
# indices: 0 1 2 3 4 0 1 2 3 4 5 6
print(" 9:", 2*w) # the list is repeated twice
w.reverse()
print("10:", w)
print("11:", w.index(7))
w.append("The last in line")
print(w)
print("12:", w.pop())
print("13:", w)
x = (1, 3, 6.62354, "Hello", "world") # example of a tuple
print(x)
print(x[3], x[4])
#x[1] = 2 # an error occurs, tuple as not mutable
y = {2,3,4,2,3, "hello", "hello"} # example of a set
print(y)
# print(y[2]) # error: 'set' object does not support indexing
z = {'a': 5, 'b': 3.1415, 'c': "some text"} # example of a dictionary
print(z)
print(z['c'])
""" Flow control """
from math import sqrt
a = 1; b = 2; c = -1;
Δ = b**2-4*a*c
print ("Δ = ", Δ)
if Δ <0:
print ("Δ is negative, no real solutions")
elif Δ==0: # note: '==' means comparison, '=' is assignment
print ("Δ is equal to 0, one real solution x =", -b/(2*a))
else:
print ("Δ is positive, two real solutions")
x1 = (-b+sqrt(Δ))/(2*a)
x2 = (-b-sqrt(Δ))/(2*a)
print ("x1 = ", x1, "x2 =", x2)
""" Loops """
for n in range(2, 15, 3):
print (n)
F = [1, 1]
for n in range(2, 10):
F.append(F[n-2]+F[n-1])
print ("Fibbonaci series: ", F)
a = 120
divisors =[]
for k in range(2,a+1):
if a%k==0:
divisors.append(k)
print ("Divisors of ", a, ": ", divisors)
""" control and some more """
for n in range(10):
if (n%2==0): continue # skip the rest and go to the beginning of the loop
if (n==7): break # finish the loop when n is equal to 7
print(n)
lst=["This", 'is', "a", 3.1415, "complicated list"]
for idx, element in enumerate(lst):
print("index=", idx, "value=", element)
squares = [k**2 for k in range(10)]
print (squares)
"""Iterations"""
x, k = 0.3, 4
for n in range(20):
x = k*x*(1-x)
print(x)
Consider a logistic map $\scriptsize x\ni[0,1]\to f(x)=kx(1-x)\in[0,1]$ for $\scriptsize k\in[0,4]$.
The iteration
$$\scriptsize x_{n+1}=kx_n(1-x_n)$$has two fixed points (solution to $\scriptsize x=f(x)$): $\scriptsize x=0$ and $\scriptsize x=1-\frac1k$.
Around a fixed point $\scriptsize x_*=f(x_*)$:
$$\scriptsize x_n=x_*+h_n$$$$\scriptsize h_{n+1}=x_{n+1}-x_*=f(x_*+h_n)-x_*$$Linear stability
$$\scriptsize h_{n+1}\approx f(x_*)+h_nf'(x_*)-x_*+\mathcal{O}(h_n^2)$$Neglecting nonlinear terms we have a geometric series:
$$\scriptsize h_{n+1}=h_nf'(x_*)$$Fixed point $x_*$ of an iteration $x_{n+1}=f(x_n)$ is linearly stable when $|f'(x_*)|<1$ and is unstable when $|f'(x_*)|>1$. Moreover the smaller $|f'(x_*)|$ the faster the convergence.
$$\scriptsize f'(0)=k\qquad f'\left(1-\frac1k\right)=2-k$$It is interesting from the point of view of a dynamical systems but very bad as a method for solving equations. But we can influence the stability
$$\scriptsize \biggl.x = f(x)\;\biggr|\;+\omega x$$$$\scriptsize (1+\omega)x = f(x)+\omega x$$$$\scriptsize x = \frac{f(x)+\omega x}{1+\omega}\equiv\tilde f(x)$$Parameter $\omega$ can control the stability ad convergence rate
$$\scriptsize \tilde f'(0)=\frac{k+\omega}{1+\omega}\qquad \tilde f'\left(1-\frac1k\right)=\frac{2-k+\omega}{1+\omega}$$"""Increase stabililty of the iteration method """
x, k, ω = 0.5, 4, 1.5
for n in range(20):
x = ( k*x*(1-x) + ω*x) / (1+ω)
print(x)
generates the following conditions (decreasing but positive sequence)
$$\scriptsize 0<I_n<I_{n-1}$$and
$$\scriptsize I_0=1-\frac{1}{e}\qquad I_n=1-nI_{n-1}$$import numpy as np
I = 1-1/np.e
for n in range(1,25):
if n>15: print("{:2d}\t{:15.8f}".format(n, I))
I = 1-n*I
""" Functions """
def Function():
print("This is a simple function. Just some code")
Function(); Function() # function calls
def Sum(a, b): # arguments
return a*b; # returned value
result = Sum(3, 222)
print(result, Sum(3,4))
def g(a, b):
return (a+b, a*b, a/b, a**b) # return a tuple of four values
s, m, d, p = g(3, 4)
print(s, m, d, p)
tup = g(3, 4)
print(tup)
Implementation of the recursive definition $$\scriptsize n!=\left\{ \begin{array}{cl} 1&\text{ for } n=0\\ n\,(n-1)!&\text{ otherwise } \end{array} \right.$$
""" Recursive functions """
def Factorial(n):
if n<1: return 1
else: return n*Factorial(n-1)
Factorial(4)
""" Recursive functions """
def Factorial(n):
print ("starting the calculation of ", n, "!")
if n<1: result = 1
else: result = n*Factorial(n-1)
print ("finished the calculation of ", n, "! = ", result)
return result
Factorial(4);
Fibonacci series: $F_1=F_2=1$ and $F_{n}=F_{n-1}+F_{n-2}$ for $n>2$
def Fib(n):
if n>2: return Fib(n-1)+Fib(n-2)
else: return 1
Fib(10)
"""Very bad Fibbonacci"""
counter = 0
def Fib(n):
global counter # global counter
counter+=1
if n<3: return 1
else: return Fib(n-1)+Fib(n-2)
print("Fib(10) = ", Fib(10))
print("counter = ", counter)
""" Functions as arguments """
import math as mt
def make_table(f):
print("==========",f.__name__,"============")
for n in range(10):
#print(f"{n*0.1:.2f}\t\t{f(n*0.1): .8f}")
print("{:.2f}\t\t{: .8f}". format(n*0.1, f(n*0.1)))
make_table(mt.sin) # another function as an argument
make_table(lambda x: x**2-1) # anonymous lambda function
"""Functions with default and key arguments """
def f1(a, b=3, c=4):
print("f1:", a, b, c)
f1(3, 4)
f1(2, c=7)
def f2(a, *b): # additional arguments are used as a tuple
print("f2:", a, b)
f2("Hello")
f2("Hello", "world", 777)
def f3(**k): # all arguments are a dictionary
print("f3:", k)
for n in k:
print("f3: ", n," = ", k[n])
f3(a=1, b=2, word="Hello")
import numpy as np # np is used as a namespace
x_list = list(range(5))
print("list: ", x_list)
x = np.arange(5)
print("np.array:", x)
print("2*x_list:", 2*x_list) # lists are joined
x = np.arange(5)
print("2*x: ", 2*x) # arrays are calculated elementwise
print("cubes: ", x**3)
x = np.linspace(0, np.pi, 8) # horizontal vector
print ("shape of x:", np.shape(x))
print ("\ntable of sins: ", np.sin(x))
y = np.array( [x, np.sin(x), np.cos(x), np.sin(x)**2+np.cos(x)**2 ] )
# 4 rows of horizontal vectors
print ("\nshape of y", np.shape(y))
print (y.transpose()) # transpose for nice output
sums=np.zeros(4)
print ("\nsums after initiation:", sums)
for k in range(4):
sums[k] = np.sum(y[:,k])
print ("Calculated sums:", sums)
Recall the example: $\scriptsize\displaystyle S = \sum_{n=1}^{N}n$.
Numpy can be much much faster than standard Python loops
%%time
N=int(1e7)
s = 0
for n in range(1,N+1):
s += n
%%time
N=int(1e7)
s = np.sum(np.arange(1,N+1))
"""Broadcasting examples"""
a = np.array([1, 2, 3])
3*a # 3*a_i
a*a # a_i^2
b = np.ones([2,3])
2*b
a*b
b*a
import matplotlib.pyplot as plt # ploting library
def make_plot():
x = np.linspace(0, 2*np.pi, 1000)
y = np.array( [ np.sin(x), np.cos(x), np.exp(-0.5*x)*np.sin(5*x) ] ).transpose()
plt.figure(figsize=(8,4), dpi=120)
plt.plot(x, y)
plt.legend(['$\sin(x)$', '$\cos(x)$', '$e^{-x/2}\sin(5x)$'])
plt.title('Example plots')
plt.xlabel('$x$')
plt.ylabel("$y$")
plt.grid(True)
plt.show()
make_plot()
"""Sympy example"""
import sympy as sp
x, k = sp.var('x k')
f = k*x*(1-x)
print("Expression:", f)
eq1 = sp.Eq(f-x)
sols = sp.solve(eq1,x)
print("Fixed points:", sols)
df = sp.diff(f,x)
print("derivative:", df)
for (n,s) in enumerate(sols):
print("f'(x_",n,") = ", sp.simplify(df.subs(x, s)))
g=f.subs(x,f) # g(x) = f(f(x))
print(g, "\n")
two_cycle = sp.solve(sp.Eq(g-x),x)
print(*two_cycle, "\n", sep="\n")
print( sp.latex(two_cycle[2]))
print(sp.latex(two_cycle[3]))