This is not an official course offered by Boston University. The sole aim of this page is to share the knowledge of how to implement Python in numerical methods. This page contains Python code for examples presented in the Fall 2015 course website.
Course Description
This course offers an advanced introduction to numerical methods for solving linear ordinary and partial differential equations, with computational implementation in Python.
Python is one of high-level programming languages that is gaining momentum in scientific computing. To work with Python, it is very recommended to use a programming environment. There are many Python's Integrated Development Environments (IDEs) available, some are commercial and others are free and open source. Personally, I would recommend the following IDEs:
Spyder: a free open-source IDE that provides MATLAB-like features, such as iPython console that works like MATLAB's command window, variable explorer which displays variables and updates statistical calculations for each variable just like MATLAB's workspace. I'd suggest installing Spyder via Anaconda.
PyCharm: also free and open source. I find PyCharm more convenient to use for working with animation and generating dynamic images, as well as for debugging a code.
All Python code available on this page were written using Python version 3.6 and higher.
To speed up Python's performance, usually for array operations, most of the code provided here use NumPy, a Python's scientific computing package. However, for comparison, code without NumPy are also presented. To see the costs of running code with different styles of coding/implementation, we compare three different ways of calculating the sum of x2 with x going from 0 to N−1 and time the execution for each method using the timeit module.
import timeit
import numpy as np
def forloopmethod(N):
a1 = [0]*N
for i in range(N):
a1[i] = float(i)**2
return sum(a1)
def listcompmethod(N):
a2 = [float(i)**2 for i in range(N)]
return sum(a2)
def numpymethod(N):
a3 = np.sum(np.arange(0, N)**2, dtype='d')
return a3
and when N=10000000, using the timeit module to time each method execution:
print(timeit.timeit('[func(N) for func in (forloopmethod,)]',
globals=globals(), number=1))
print(timeit.timeit('[func(N) for func in (listcompmethod,)]',
globals=globals(), number=1))
print(timeit.timeit('[func(N) for func in (numpymethod,)]',
globals=globals(), number=1))
we obtain that the first method with conventional for loop takes the longest time to run, about 3.4 seconds, the second method with list comprehension runs for about 3.2 seconds, and the method with NumPy runs the fastest, about 0.1 seconds.
An example code to measure execution time is available here.
Lecture 7:
This lecture discusses different numerical methods to solve ordinary differential equations, such as forward Euler, backward Euler, and central difference methods. Below are simple examples of how to implement these methods in Python, based on formulas given in the lecture note (see lecture 7 on Numerical Differentiation above). We also learn how to pass multiple arguments using the magic variable with the asterisk (*) symbol.
In this extra handout for lecture 8 [pdf], details on how to create functions in Python for the following basic Euler methods are discussed. These methods are used to solve:dydx=3(1+x)−yanddydx=2x−4xy
The implementation of Runge-Kutta methods in Python is similar to the Heun's and midpoint methods explained in lecture 8. Here we discuss 2nd-order Runge-Kutta methods with A=12 (type A), A=0 (type B), A=13 (type C), as well as 3rd-order, 4th-order, and Runge-Kutta-Fehlberg (RKF45) methods. These methods are used to solve the following ODE,dydx=e−2x−2y
with initial condition y(0)=0.1, interval x=[0,2], and step size h=0.2.
In the files above, these methods are used to solve:dydx=x−y2with x=[0,3], y(0)=1.0, and h=0.125. These methods need to invoke other methods, such as Runge-Kutta methods, to get their initial values.
Systems of ODEs, such as the Van der Pol oscillatordy1dt=y2anddy2dt=a(1−y21)y2−y1and a system of activator-inhibitor equationsdxdt=a+bx21+x2+ry−xanddydt=ε(cx+y0−y)where both result in oscillating solutions, need to be solved with high accuracy solvers. We use the following methods:
Python code for these methods from previous lectures can be directly used for multiple ODEs, except for the 4-step Adams-Bashforth-Moulton method, where we need to modify the variable yn = yy[0:m] and several variables within the for loop (highlighted in blue):
We employ a second-order finite difference formula to solve the following boundary value problem (BVP):d2ydx2=12x2for x=[0,1] with y(0)=0 and y(1)=0. The exact solution of this problem is y(x)=x4−4.
This lecture discusses how to numerically solve the Poisson equation,−∇2u=fwith different boundary conditions (Dirichlet and von Neumann conditions), using the 2nd-order central difference method. In particular, we implement Python to solve,−∇2u=20cos(3πx)sin(2πy)
Dirichlet problem: poissonDirichlet.py with boundary conditions u(0,y)=y2,u(1,y)=1,u(x,0)=x3,u(x,1)=1.
Von Neumann problem: poissonNeumann.py with boundary conditions ux(0,y)=0,ux(1,y)=0,uy(x,0)=0,uy(x,1)=0.
This lecture discusses how to numerically solve the 1-dimensional reaction-diffusion equation,∂u∂t=D∂2u∂x2+αuusing forward time central space (FTCS), backward time central space (BTCS), and Crank-Nicolson (CN) methods.
Periodic boundary conditions are applied in all methods.
*All code have been modified for animation/dynamic simulations, by removing plt.show() from the class loop and placing it inside of the main() function.
This lecture discusses how to numerically solve the 2-dimensional diffusion equation,∂u∂t=D∇2uwith zero-flux boundary condition using the ADI (Alternating-Direction Implicit) method.
Fig. Solution of 2D diffusion equation using the ADI method.
Lecture 21:
In this lecture, we solve the 2-dimensional wave equation,∂2u∂t2=D(∂2u∂x2+∂2u∂y2)using:
The finite difference method, by applying the three-point central difference approximation for the time and space discretization. This way of approximation leads to an explicit central difference method, where it requiresr=4DΔt2Δx2+Δy2<1to guarantee stability. The following example is a solution of the wave equation on a [0,2]×[0,2] domain, with diffusion coefficient D=0.25, initial condition u(x,y,0)=0.1sin(πx)sin(πy2), initial velocity ∂u(x,y,0)∂t=0, and Dirichlet boundary condition u(0,y,t)=u(2,y,t)=u(x,0,t)=u(x,2,t)=0.
For the requirement of r<1, we use Python's assert statement, so that the program will not execute and raise an error if the requirement is not fulfilled. In the code below, the assertion is applied in the initialization function. Try running the code with higher diffusion coefficient, such as D=1.5, by modifying simulator = WaveEquationFD(200, 0.25, 50, 50) to simulator = WaveEquationFD(200, 1.5, 50, 50) and see how the assertion works.
FD_wave1.py Fig. Solution of 2D wave equation using finite difference method.
A Spectral method, by applying a leapfrog method for time discretization and a Chebyshev spectral method on a tensor product grid for spatial discretization. The following example is a solution of the wave equation on a [−1,1]×[−1,1] domain, with diffusion coefficient D=1.0, initial condition u(x,y,0)=exp(−40((x−0.4)2+y2)), initial velocity ∂u(x,y,0)∂t=0, and Dirichlet boundary condition u(−1,y,t)=u(1,y,t)=u(x,−1,t)=u(x,1,t)=0. This method uses a computational spectral grid, clustered at the boundaries. Chebyshev differentiation is carried out by the fast Fourier transform. The results at each grid point are spectrally accurate, despite errors of magnitude O((Δt)2) caused by time-stepping[1].
Spectral_wave2.py Fig. Solution of 2D wave equation using a spectral method.
[1] Lloyd N. Trefethen, Spectral Methods in MATLAB, SIAM.
No hay comentarios:
Publicar un comentario