• What is Python?
• Some general facts about Python
• Programming 101
• Selected packages for scientific computation:
  NumPy, SciPy, matplotlib, SymPy, FEniCS/DOLFIN

• free software
• universal programming language
• many packages existing
• easily extendable

Python is (line many script languages) very intuitive.

• no explicit type specifications (duck typing)
• code blocks must be indented (leads to improved readability)
• Very high-quality code.
  0.005 defects per thousand lines of code:
  "The open-source Python programming language [...] now surpasses that of its open-source and proprietary peers, according to a study published by development testing vendor Coverity." -- Sean M. Kerner, eweek.com

• Holistic Language Design
  "Doing mathematical programming in a general-purpose language is easier than doing general-purpose programming in a mathematical language."
• Readability
• Balance of High Level and Low Level Programming
• Language Interoperability
• Documentation System
• Hierarchical Module System
• Data Structures
• Available Libraries
• Testing Framework
• Large user base

• Not as popular as, e.g., MATLAB(r)
  
• ...
  

print('Hello world!')
a = 42                           # integer
b = 42 / 21 + 10                 # = 12
my_list = [1, 3, 'Hello World!'] # list: all types permitted!
print(my_list[0])                # 0-based indexing
for i in range(4):               # i in 0,...,3
    print(i)
for elem in my_list:
    print(elem)
d = { 'Hallo': 'hello',       # dictionary ("map")
      'Welt': 'world',
      42: 'fourty-two'}
print(d['Hallo'])             # -> hello
tuple = (4,3,1)
if tuple == (4,):             # easy to compare!
    print('yeah')

# mymodule.py
def square(a):       # functions can be defined anywhere
    return a*a       # with arbitrary names (like variables)

def arnoldi(A, v, maxiter=10, ortho='mgs'):  # default arguments
    # compute V, H
    return V, H      # multiple return arguments

import mymodule      # import is possible anywhere in the code
print( mymodule.square(3) )

from mymodule import square, arnoldi
print( square(3) )
V, H = arnoldi(A, v, ortho='house') # maxiter is set to default (10)

NumPy is the basis package for numerics (based on BLAS, LAPACK, etc.).

from numpy import array, ones, dot   # N-dim Arrays
v = array([1,2,3])           # 1-dim arr, shape==(3,)
v = ones(5)                  # 1-dim arr, shape==(5,)
v = ones((5,1))              # 2-dim arr, shape==(5,1)
A = array([[1,2,3],[4,5,6]]) # 2x3-matrix: 2-dim arr, shape==(2,3)
A = ones((5,5))              # 5x5-matrix: 2-dim arr, shape==(5,5)
b = dot(A, v)                # A*v

The default data type in NumPy are NumPy arrays; there are also NumPy matrices. The main difference is that one can use the *-operator for matrix-matrix multiplication (as in MATLAB). Many (especially the NumPy community) recommend NumPy arrays.

Well-known commands in NumPy:

• zeros, ones, diag, abs, linspace, ...
• module random: rand, randn, ...
• module linalg: solve, norm, eig, svd, lu, qr, cholesky, ...
• module polynomial: monomial, Chebyshev, Legendre, ...

MATLAB's \-operator: solve(A,b).

For more info, see NumPy documentation.

Interesting for us:

• scipy.linalg: extended functionality: banded matrices, qz, matrix functions, ...
• scipy.sparse: several sparse matrix types, cg, gmres, eigs, ...
• scipy.optimize: simplex, newton, leastsq, ...

For more info, see SciPy documentation.

from matplotlib import pyplot as pp

data = [x**2 for x in range(100)] # list comprehension!
pp.plot(data)
pp.show()

LaTeX-integration via matplotlib2tikz.

import sympy as smp

x, y = smp.symbols('x, y')

f = x**2 + sin(x) * y
diff(f, x)

diff(f, x, 2)

• differentiation
• integration
• ...

• Powerful interactive shells (terminal and Qt-based).
• A browser-based notebook with support for code, text, mathematical expressions, inline plots and other rich media.
• Support for interactive data visualization and use of GUI toolkits.
• Flexible, embeddable interpreters to load into your own projects.
• Easy to use, high-performance tools for parallel computing.

• krypy (by André Gaul, TU Berlin) - Krylov subspace methods & tools

• pyamg - algebraic multigrid

• ...

J = [ X(NODES(2),1)-X(NODES(1),1) , X(NODES(3),1)-X(NODES(1),1)  ;
      X(NODES(2),2)-X(NODES(1),2) , X(NODES(3),2)-X(NODES(1),2)  ];
AbsDetJ = abs( J(1,1)*J(2,2) - J(1,2)*J(2,1) ) ;
L = (eps/(2*AbsDetJ)) * [ J(1,2)^2 + J(2,2)^2  , - J(1,1)*J(1,2) - J(2,1)*J(2,2)  ;
                                     0         ,   J(1,1)^2 + J(2,1)^2       ];
A = z3A(1,1) =   L(1,1) + 2*L(1,2) + L(2,2) + AbsDetJ/12 ;
A(1,2) = - L(1,1) - L(1,2)            + AbsDetJ/24 ;
A(1,3) = - L(1,2) - L(2,2)            + AbsDetJ/24 ;

A(2,1) = - L(1,1) - L(1,2)            + AbsDetJ/24 ;
A(2,2) =   L(1,1)                     + AbsDetJ/12 ;
A(2,3) =   L(1,2)                     + AbsDetJ/24 ;

A(3,1) = - L(1,2) - L(2,2)            + AbsDetJ/24 ;
A(3,2) =   L(1,2)                     + AbsDetJ/12 ;
A(3,3) =   L(2,2)                     + AbsDetJ/12 ;
    

for r=1:size(ECT,1)
  localStiffnessMatrix = local_assemble(ECT(r,:),X,eps);
  for alpha=1:3
     i = ECT(r,alpha);
     for beta= 1:3
        j = ECT(r,beta);
        K(i,j) = K(i,j) + localStiffnessMatrix(alpha,beta);
     end
  end
end
    

• Quadrature of RHS
• boundary conditions
• ...

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'Lagrange', 1)

# Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')
def u0_boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, u0, u0_boundary)
    

[...]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution and mesh
plot(u)
plot(mesh)

# Dump solution to file in VTK format
file = File('poisson.pvd')
file << u

# Hold plot
interactive()
    

prm = parameters['krylov_solver'] # short form
prm['absolute_tolerance'] = 0.0
prm['relative_tolerance'] = 1E-13
prm['maximum_iterations'] = 10000
prm['monitor_convergence'] = True
solve(a == L, u, [bc0, bc1],
      solver_parameters={'linear_solver': 'cg',
                         'preconditioner': 'ilu'}
      )
    

    from dolfin import *
    parameters['linear_algebra_backend'] = 'uBLAS'
    [...]
    a = inner(nabla_grad(u), nabla_grad(v))*dx
    L = f*v*dx + g*v*ds

    A = assemble(a, bcs=[bc0,bc1])
    b = assemble(L, bcs=[bc0,bc1])
    

$ python ex5.py
(array([    0,     3,     8, ..., 11433, 11438, 11441]), array([   0,    1,    2, ..., 1678, 1679, 1680]), array([ 1. , -0.5, -0.5, ...,  0. ,  0. ,  1. ]))
[ 0.00114583  0.0065625   0.         ...,  0.          1.          1.        ]
    

• Build your mesh: Gmsh, tetgen, CUBIT, ... (formats: VTK, Exodus,...)
• Feed the mesh into the solver.

UnitSquare(4, 4)
UnitCube(4, 4, 4)
Rectangle(a, 0, b, 1, nr, nt, 'crossed')
[...]
    

    $ dolfin-convert mymesh.msh mymesh.xml
    dolfin-convert mymesh.msh mymesh.xml
    Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format
    Expecting 2001 vertices
    Found all vertices
    Expecting 10469 cells
    Found all cells
    Conversion done

    $ ls mymesh*.xml
    mymesh_physical_region.xml  mymesh.xml
    

    a = u*v*dx + dt*inner(grad(u), grad(v))*dx
    L = y_1*v*dx + dt*f*v*dx
    while t < T:
        solve(a == L, y)
        t += dt
        y_1.assign(y)