SymPy - Python library for symbolic mathematics

http://sympy.org/en/index.html

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python and does not require any external libraries.

matplotlib: python plotting

http://matplotlib.sourceforge.net/

matplotlib is a python 2D plotting library which produces publication quality figures in a variety of hardcopy formats and interactive environments across platforms. matplotlib can be used in python scripts, the python and ipython shell (ala MATLAB®* or Mathematica®†), web application servers, and six graphical user interface toolkits. matplotlib tries to make easy things easy and hard things possible. You can generate plots, histograms, power spectra, bar charts, errorcharts, scatterplots, etc, with just a few lines of code. For a sampling, see the screenshots, thumbnail gallery, and examples directory

SciPy and NumPy - Scientific Tools for Python

Scientific Tools for Python 

SciPy (pronounced "Sigh Pie") is open-source software for mathematics, science, and engineering. It is also the name of a very popular conference on scientific programming with Python. The SciPy library depends on NumPy, which provides convenient and fast N-dimensional array manipulation. The SciPy library is built to work with NumPy arrays, and provides many user-friendly and efficient numerical routines such as routines for numerical integration and optimization. Together, they run on all popular operating systems, are quick to install, and are free of charge. NumPy and SciPy are easy to use, but powerful enough to be depended upon by some of the world's leading scientists and engineers. If you need to manipulate numbers on a computer and display or publish the results, give SciPy a try!

Python and Scientific Computing 

NumPy and SciPy are two of many open-source packages for scientific computing that use the Python programming language. This website, together with other subdomains of the scipy.org domain, serves as a portal for all scientific computing with Python, not just NumPy and SciPy. The index under Topical Software in the navigation bar lists these domains and other destinations for scientific software using Python.  

Good places to start to learn more about SciPy: 

• More about NumPy / SciPy

• scikits add-on toolkits for numpy/scipy

• Index of Python packages for scientific computing

• Projects using NumPy / SciPy

• Download NumPy / SciPy

• Install SciPy

• Documentation

• Cookbook

• Frequently Asked Questions

This is a community effort. We seek volunteers at all levels of ability to work on the project, from coding and packaging to documentation, tutorials, recipes, and the web site. Visit the Developer Zone if you are interested in helping out. SciPy is sponsored by Enthought, Inc.

SciPy optimization on Ubuntu Linux

NumPy is the fundamental package needed for scientific computing with Python. It contains among other things: 

• a powerful N-dimensional array object 

• sophisticated (broadcasting) functions 

• tools for integrating C/C++ and Fortran code 

• useful linear algebra, Fourier transform, and random number capabilities. 

Besides its obvious scientific uses, NumPy can also be used as an efficient multi-dimensional container of generic data. Arbitrary data-types can be defined. This allows NumPy to seamlessly and speedily integrate with a wide variety of databases. Numpy is licensed under the BSD license, enabling reuse with few restrictions.

OpenFOAM and installation of ver. 2.1.0 in Ubuntu 11.10

I've been wanting to install and try this software for awhile now.  It is called OpenFOAM and is a free and open source CFD and numerical software.  While ANSYS Fluent and others are great and heavily developed, they are quite expensive and un-readily available to many users.  I've done some research and OpenFOAM seems to be the leader of the free world.  They seem to be very open and constantly updated and developed.

About the OpenFOAM Foundation: http://www.openfoam.org/index.php

The OpenFOAM® Foundation is a non-stock corporation, incorporated on 5th August 2011 in Delaware, USA. It is a nonprofit organization whose “specific objectives and purposes... shall be to promote and manage the free open source distribution of the OpenFOAM software” (Bylaws, Article II). The Foundation has taken on the guardianship of the OpenFOAM software, previously undertaken by OpenCFD Ltd. Through its bylaws, it formalises the commitment, begun by OpenCFD in 2004, to ensure OpenFOAM is free and open source only. It is supported by SGI and lists Mark Barrenechea (CEO, SGI) and Henry Weller (creator of OpenFOAM) amongst its Board of Directors. 

The Foundation distributes the current repository version of OpenFOAM (and will distribute future ‘version’ releases of OpenFOAM) under the GNU general public licence (GPL). The GPL gives users the freedom to modify and redistribute the software and a guarantee of continued free use — as long as the terms of the GPL are adhered to. There are two elements to the GPL that provide protection against exploitation by companies including OpenFOAM within non-free and/or closed source software products. First, when any modified version of the software is redistributed, the source code must also be made available by the distributor. Secondly, any modified version can only legally be distributed open source under the GPL and software that links intimately enough to OpenFOAM has to be distributed under the GPL as well. 

People and organisations contribute to OpenFOAM on the understanding that their contributions become part of the public commons of free software and are not exploited by producers of non-free software. The Foundation exists to ensure this, adopting a suitably strong, open source license, i.e. the GPL. It will manage the software code base and promote the software, e.g. through its website www.openfoam.org, whose content will grow over time. OpenFOAM Foundation 15th August 2011

About OpenFOAM: http://www.openfoam.com/

The OpenFOAM® (Open Field Operation and Manipulation) CFD Toolbox is a free, open source CFD software package produced by OpenCFD Ltd. It has a large user base across most areas of engineering and science, from both commercial and academic organisations. OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to solid dynamics and electromagnetics. It includes tools for meshing, notably snappyHexMesh, a parallelised mesher for complex CAD geometries, and for pre- and post-processing. Almost everything (including meshing, and pre- and post-processing) runs in parallel as standard, enabling users to take full advantage of computer hardware at their disposal. 

By being open, OpenFOAM offers users complete freedom to customise and extend its existing functionality, either by themselves or through support from OpenCFD. It follows a highly modular code design in which collections of functionality (e.g. numerical methods, meshing, physical models, …) are each compiled into their own shared library. Executable applications are then created that are simply linked to the library functionality. OpenFOAM includes over 80 solver applications that simulate specific problems in engineering mechanics and over 170 utility applications that perform pre- and post-processing tasks, e.g. meshing, data visualisation, etc.

The OpenFOAM Team: http://www.openfoam.com/about/

The OpenFOAM Team produces the OpenFOAM® open source CFD toolbox and distributes it through the OpenFOAM Foundation. The OpenFOAM Team was formed at OpenCFD Ltd, a company established in 2004 to coincide with the release of its OpenFOAM software under general public licence. The company was founded by: Henry Weller, the creator of OpenFOAM and its chief developer from inception to the present day; Chris Greenshields; and, Mattijs Janssens. In 2011, the team moved to SGI Corp, following its acquisition of OpenCFD Ltd.

OpenFOAM Services:

The OpenFOAM Team has unrivalled experience with OpenFOAM. Over the past 7 years they have successfully provided services to numerous science/engineering companies, consultancies and universities. Developments: written over half a million of lines of code. Support: delivered over 10,000 hours of support. Training: delivered over 100 training courses to over 1000 OpenFOAM users. The OpenFOAM Team is Henry Weller, Chris Greenshields, Mattijs Janssens, Andy Heather, Sergio Ferraris, Jenya Collings, Gijs Wierink and Laurence McGlashan.

I recommend you check this out if you are interested in a free CFD simulator.  While it may not have all the bells and whistles of Fluent or other similar software, but I bet you can get just as much done and in the process understand what you are doing much better as it can be the tendency to treat more fancy software as a black box.  OpenFOAM has meshing and post processing capabilities (ParaView), supports many imports of meshes from other software, has many solvers, many models, and many libraries, it also has great documentation, tutorials, and community support.  Training and customer support is also offered.

To install OpenFOAM simply go to http://www.openfoam.org/download/ where they have choices of the Ubuntu Deb Pack, SuSE RPM Pack, the Source Pack, and the Git Repository.  Luckily for us, I am using Ubuntu 11.10 so I simply went to their page for that installation here: http://www.openfoam.org/download/ubuntu.php

I followed the instructions and had no problems until I got to the end where it says Getting Started.  It seems to assume that a directory exist in your home directory called OpenFOAM.  However, this is not where OpenFOAM was installed as you can see it is installed in /opt/openfoam210 instead.  So you need to first create the directory OpenFOAM in your home folder then everything else should be fine in the directions.

This is the link for the documentation (tutorials and more): http://www.openfoam.org/docs/user/

One downside to OpenFOAM is that it is not cross platform.  It only works on Linux such as Ubuntu, Fedora, SuSE, etc.

Here is a screen shot showing that OpenFOAM installed into the opt directory.

Here is a screen shot showing the /home/timbarber/OpenFoam/timbarber-2.1.0/run/tutorial tree.  This is the step where the tutorial material is copied from the /opt/openfoam210 directory.

Here is the message I get after editing the .bashrc file then either typing source $HOME/.bashrc and then typing icoFoam -help or opening a new terminal window then typing icoFoam -help

Screen shot after running blockMesh in the cavity directory.

Viewing the mesh in ParaView by typing paraFoam in the terminal.

After following the instructions on OpenFOAM's user document web page, this is a view of the mesh.

Numerical Math - Solving systems of equations utilizing matrices (data fitting)

A nice pointer from Laundau et al [1] in the opening of their chapter (Ch. 8) on Solving Systems of Equations with Matrices; Data Fitting.  Many physical systems are modeled utilizing matrices which consist of a system of simultaneous equations.  However, these sets can become quite large and complicated which is why computers are very good with these processes.  Usually the algorithms for solving these sets of equations that utilize matrix theory is that they require repetition of a small list of steps which have been written in an efficient method.

One additional technique for speed is to tune the algorithm to the actual architecture of the computer which Landau et al [1] discuss more in their Chapter 14 called High-Performance Computing Hardware, Tuning, & Parallel Computing.

Many libraries exist which are "industrial-strength" subroutines for solving these matrix systems.  A majority of these libraries are well established such as the IMSL Numerical Libraries by Rogue Wave Software, Inc., the GNU Scientific Library (GSL), and the Netlib Repository at UTK and ORNL which contains LAPACK — Linear Algebra PACKage.  Landau et al [1] note that these libraries are usually an order of magnitude faster than general methods in linear algebra texts.  The libraries are streamlined for minimal round-off error and are aimed to solve a large spectrum of problems with high success.  It is here and for the reasons just mentioned where Landau et al warn that it is best if you don't write your own matrix subroutines but retrieve them from one of these libraries.  The libraries also provide the advantage of allowing the user to run them on one machine/processor or many machines/processors by varying with the computer architecture.

Next, Landau et al [1] ask the question which proposes to the user what is considered a "large" matrix.  Before, a large matrix was based upon a fraction of the RAM available to the computer system.  However, Landau et al describe a "large" matrix as now based upon the numerical time it takes to obtain values.  That is, if any waiting time is required, then a library should be used.  Landau et al also comment that the libraries are beneficial for speed even when the matrices may be small (which might apply in graphics processing).

One negative side effect lies in the multiple languages that the libraries are written.  That is one library may be in Fortran while the user is a C coder.  However, today libraries might exist which are programmed in or for many different computer coding languages.

---

References:


[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

deal.II Homepage

deal.II Homepage: A Finite Element Differential Equations Analysis Library - October 9th, 2011: deal.II 7.1 released

What is deal.II? 

deal.II is a C++ program library targeted at the computational solution of partial differential equations using adaptive finite elements. It uses state-of-the-art programming techniques to offer you a modern interface to the complex data structures and algorithms required. 

The main aim of deal.II is to enable rapid development of modern finite element codes, using among other aspects adaptive meshes and a wide array of tools classes often used in finite element program. Writing such programs is a non-trivial task, and successful programs tend to become very large and complex. We believe that this is best done using a program library that takes care of the details of grid handling and refinement, handling of degrees of freedom, input of meshes and output of results in graphics formats, and the like. Likewise, support for several space dimensions at once is included in a way such that programs can be written independent of the space dimension without unreasonable penalties on run-time and memory consumption. 

deal.II is widely used in many academic and commercial projects. For its creation, its principal authors have received the 2007 J. H. Wilkinson Prize for Numerical Software. It is also part of the industry standard SPEC CPU 2006 benchmark suite used to determine the speed of computers and compilers, and comes pre-installed on the machines offered by the commercial Sun Grid program. 

deal.II emerged from work at the Numerical Methods Group at Universität Heidelberg, Germany, which is at the forefront of adaptive finite element methods and error estimators. Today, it is maintained by two of its original authors at Texas A&M University, and dozens of contributors and several hundred users are scattered around the world (see our credits page for a detailed list of people contributing to deal.II).

What deal.II can offer you? 

If you are active in the field of adaptive finite element methods, deal.II might be the right library for your projects. Among other features, it offers: 

Support for one, two, and three space dimensions, using a unified interface that allows to write programs almost dimension independent. 

Handling of locally refined grids, including different adaptive refinement strategies based on local error indicators and error estimators. Both h, p, and hp refinement is fully supported for continuous and discontinuous elements. 

Support for a variety of finite elements: Lagrange elements of any order, continuous and discontinuous; Nedelec and Raviart-Thomas elements of any order; elements composed of other elements. 

Parallelization on single machine through the Threading Build Blocks and across nodes via MPI. deal.II has been shown to scale to at least 16k processors. 

Extensive documentation: all documentation is available online in a logical tree structure to allow fast access to the information you need. If printed it comprises more than 500 pages of tutorials, several reports, and presently some 5,000 pages of programming interface documentation with explanations of all classes, functions, and variables. All documentation comes with the library and is available online locally on your computer after installation. 

Modern software techniques that make access to the complex data structures and algorithms as transparent as possible. The use of object oriented programming allows for program structures similar to the structures in mathematical analysis. 

A complete stand-alone linear algebra library including sparse matrices, vectors, Krylov subspace solvers, support for blocked systems, and interface to other packages such as Trilinos, PETSc and METIS.

Support for several output formats, including many common formats for visualization of scientific data.

Portable support for a variety of computer platforms and compilers.
Free source code under an Open Source license, and the invitation to contribute to further development of the library.

Numerical Math - The Jacobi method

A simple explanation giving by Landau et al [1] defines the Jacobi method as an initial sweep of a values without updating if available. The Jacobi method is very basic. For example, from the square wire problem in Landau et al the numerical algorithm to be solved is

\[ U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \]

The initialization and BCs symmetry are preserved in this way.

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

Numerical Math - The Gauss-Seidel (GS) method

The Gauss-Seidel method is an improvement upon the basic Jacobi method [1]. The GS method utilizes known values and current updates them in the algorithm as opposed to the Jacobi method which sweeps across the domain without utilizing known values to accelerate convergence. The accelerated convergence also produces less round-off error and utilizes less memory since only one generation of guesses is needed.

If the sweep begins at the top left corner then the GS algorithm looks like

\begin{equation} U_{i, j} = \dfrac{1}{4} \left[ U^{(old)}_{i + 1, j} + U^{(new)}_{i - 1, j} + U^{(old)}_{i, j + 1} + U^{(new)}_{i, j - 1} \right] \end{equation}

Compare to the Jacobi method.  Both algorithms shown are for the square wire finite difference technique problem in Landau et al [1].

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \end{equation}

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

Numerical Math - The Finite Element Method (FEM)

From Landau [1], the finite element method (FEM) is explained as solving a PDE where the whole region or domain is subsectioned into smaller areas known as elements. Next an initial solution is formulated for the PDE in each of the elements.  The next step comes in the form of modifying the parameters of the initial formalization known as a best fit.

Pepper and Heinrich [2] open their FEM book by explaining that the method is a numerical formulation that solves physics and engineering models described by differential equations  (DEs). Similar to finite difference methods (FDM), FEM models a geometric region which is then broken up into a number of smaller subregions which results in a network known as a mesh. One difference between FDM and FEM includes mesh type (geometrically). The FDM requires that the network be made of orthogonal rows and columns (squares and rectangles) while the FEM does not require this limitation and can be, in fact, any shape such as triangles and/or quadrilaterals in two dimensions and tetrahedrons and/or hexahedrons in three dimensions.  Next, each FE is initialized with an approximate functions of the unknown variables for which to be solved. Expansions determine the variables approximations and appear as a linear or higher-order polynomials functions. These expansions, in turn, depend upon the geometric shape of the elements and location known as nodes. A second difference between FDM and FEM noted by Pepper and Heinrich [3] entails the solution method. The FEM integrates over each subregion then adds or connects them to make up the whole. This integration and summation engenders a set of finite linear equations per section which can then be solved utilizing linear algebra methods. Jiang [6] also discusses that FEM does not...

The network of elements and nodes or where these elements connect make up discrete systems such as a trusses, circuits, and fluid transport pipes [5]. In order to solve for the system variables such as displacements, electric potentials, and pressures, one can begin with a known and simple parameter such as

Hooke’s Law:

\begin{equation} F = \dfrac{\Delta L}{R} = E A \dfrac{\Delta L}{L} \end{equation}

Ohm’s Law:

\begin{equation} i = \dfrac{\Delta V}{R} = \dfrac{A}{\rho} \dfrac{\Delta V}{L} \end{equation}

Poiseuille’s Law:

\begin{equation} \dot{m} = \dfrac{\Delta p}{R} = \dfrac{\rho \pi D^4}{128 \mu} \dfrac{\Delta p}{L} \end{equation}

Jiang [6] notes that FEM has been utilized as one of the most general numerical techniques to solve DEs and done so with great success. Jiang even quotes Oden as

“Perhaps no other family of approximation methods has had a greater impact on the theory and application of numerical methods during the twentieth century”

In progress...to be continued.

---

References:

[1] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. The Finite Element Method - It’s Basis and Fundamentals, 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

[2] D. W. Pepper and J. C. Heinrich. The Finite Element Method: Basic Concepts and Applications, Taylor & Francis Hemisphere Publishing Corporation,. 1992

[4] J. C. Heinrich and D. W. Pepper. The Intermediate Finite Element Method: Fluid Flow and Heat Transfer Applications, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1999

[5] G. Comini, S. D. Giudice and C. Nonino. Finite Element Analysis in Heat Transfer: Basic Formulation & Linear Problems, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1994

[6] B.-N. Jiang. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Springer-Verlag, Berlin, Germany. 1998

[7] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. The Finite Element Method for Fluid Dynamics, (Volume 3) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

Some references for numerical methods and CFD

A list of numerical analysis and CFD book references.

[1] O. C. Zienkiewicz, R. L. Taylor, and J. Z. Zhu. The Finite Element Method - It’s Basis and Fundamentals, (Volume 1) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

[2] O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method for Solid and Structural Mechanics, (Volume 2) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005

[3] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. The Finite Element Method for Fluid Dynamics, (Volume 3) 6th ed. Elsevier Butterworth-Heinemann, Burlington, MA. 2005. 2005. 2005. 2005

[3] D. W. Pepper and J. C. Heinrich. The Finite Element Method: Basic Concepts and Applications, Taylor & Francis Hemisphere Publishing Corporation,. 1992

[4] G. Comini, S. D. Giudice and C. Nonino. Finite Element Analysis in Heat Transfer: Basic Formulation & Linear Problems, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1994

[5] J. C. Heinrich and D. W. Pepper. The Intermediate Finite Element Method: Fluid Flow and Heat Transfer Applications, Taylor & Francis Hemisphere Publishing Corporation, Washington, DC. 1999

[6] J. Donéa and A. Huertu. Finite Element Methods for Flow Problems, John Wiley & Sons, Ltd., Chinchester, UK. 2003

[7] B.-N. Jiang. The Least-Squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics, Springer-Verlag, Berlin, Germany. 1998

[8] P. M. Gresho and R. L. Sani. Incompressible Flow and the Finite Element Method, Incompressible Flow & the Finite Element Method - Advection-Diffusion & Isothermal Laminar Flow, John Wiley & Sons, Ltd., Chinchester, UK. 1998

[9] P. M. Gresho and R. L. Sani. Incompressible Flow and the Finite Element Method, Volume 1, Advection-Diffusion, (reprint), John Wiley & Sons, Ltd., Chinchester, UK. 2000

[10] P. M. Gresho and R. L. Sani. Incompressible Flow and the Finite Element Method, Volume 2, Isothermal Laminar Flow, (reprint), John Wiley & Sons, Ltd., Chinchester, UK. 2000

[11] R. H. Landau and M. J. Páez. Computational Physics: Problem Solving with Computers, John Wiley & Sons, New York, NY. 1997

[12] R. H. Landau and M. J. Páez. Computational Physics: Problem Solving with Computers, 2nd ed. John Wiley & Sons, New York, NY. 2007

[13] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics - Introductory Computational Science, Princeton University Press, Princeton, New Jersey. 2008

[14] T. J. Chung. Computational Fluid Dynamics, Cambridge University Press,Cambridge, UK. 2002

[15] T. J. Chung. Computational Fluid Dynamics, 2nd ed. Cambridge University Press, Cambridge, UK. 2010

[16] D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. Computational Fluid Mechanics and Heat Transfer, Taylor & Francis Hemisphere, . 1984

[17] D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. Computational Fluid Mechanics and Heat Transfer, 2nd ed. Taylor & Francis  Hemisphere, New York, NY. 1997

[18] S. V. Patankar. Numerical Heat Transfer and Fluid Flow, Taylor & Francis Hemisphere, . 1980

[18] T.-M. Shih. Numerical Heat Transfer, Taylor & Francis Hemisphere, . 1984

[19] J. N. Reddy and D. K. Gartling. The Finite Element Method in Heat Transfer and Fluid Dynamics, 2nd ed. CRC Press, . 2001

[20] J. N. Reddy and D. K. Gartling. The Finite Element Method in Heat Transfer and Fluid Dynamics, 3rd ed. CRC Press, . 2010

[21] A. W. Date. Introduction to Computational Fluid Dynamics, Cambridge University Press, Cambridge, UK. 2005

[22] J. Tu, G. H. Yeoh, and C. Liu. Computational Fluid Dynamics: A Practical Approach, Elsevier Butterworth-Heinemann, Burlington, MA. 2008

[23] J. Blaz̆ek. Computational Fluid Dynamics: Principles and Applications, (reprinted in 2006) Elsevier, Oxford, UK. 2001

[24] J. Blaz̆ek. Computational Fluid Dynamics: Principles and Applications, 2nd ed. (reprinted in 2007) Elsevier, Oxford, UK. 2005

[25] H. K. Versteeg and W. Malalasekera. An Introduction to Computational Fluid Dynamics: The Finite Volume Method Approach, Longman Scientific & Technical, Harlow, UK. 1995

[26] H. K. Versteeg and W. Malalasekera. An Introduction to Computational Fluid Dynamics: The Finite Volume Method Approach, 2nd ed. Prentice Hall Pearson Education Limited, Harlow, UK. 2007

[27] J. D. Hoffman. Numerical Methods for Engineers and Scientists. 2nd ed. Marcel Dekker, Inc., New York, NY. 2001

In progress...to be continued.

Numerical Math - Square wire problem - finite difference method

The square-wire problem is described by the 2-D Laplace and Poisson equations [1].

\begin{equation} \dfrac{\partial^2 U \left( x, y \right)}{\partial x^2} + \dfrac{\partial^2 U \left( x, y \right)}{\partial y^2} = \left\{ \begin{array}{cr} 0 & \text{Laplace's Equation}\\ -4 \pi \rho \left( \mathbf{x} \right) & \text{Poisson's Equation} \end{array} \right. \end{equation}

Next, in order to obtain the second order derivatives a Taylor series expansion with steps of positive and negative values in both \( x \) and \( y \) are executed for \( U \).

\begin{equation} U \left( x + \Delta x, y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x - \Delta x, y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial x} \Delta x + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 - \cdots, \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) = U \left( x, y \right) + \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + \cdots, \end{equation}

\begin{equation} U \left( x, y - \Delta y \right) = U \left( x, y \right) - \dfrac{\partial U}{\partial y} \Delta y + \dfrac{1}{2} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 - \cdots, \end{equation}

Now, adding both equations each for \( x \) and \( y \) the odd terms cancel, and truncating at \( O \left( \Delta^4 \right) \) obtains

\begin{equation} U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 + O \left( \Delta x^4 \right) \end{equation}

\begin{equation} U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) = 2 U \left( x, y \right) + \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 + O \left( \Delta y^4 \right) \end{equation}

Rearranging to segregate the second derivative engenders

\begin{equation}  \dfrac{\partial^2 U}{\partial x^2} \left( \Delta x \right)^2 = U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right) \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} \left( \Delta y \right)^2 = U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right) \end{equation}

Finally, dividing produces

\begin{equation} \dfrac{\partial^2 U}{\partial x^2} = \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \end{equation}

\begin{equation} \dfrac{\partial^2 U}{\partial y^2} = \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} \end{equation}

Next, substituting the second order derivatives into the Poisson equation produces

\begin{equation} \begin{array}{l} \dfrac{U \left( x + \Delta x, y \right) + U \left( x - \Delta x, y \right) - 2 U \left( x, y \right)}{\Delta x^2} \\ \hspace{1.15 in} + \dfrac{U \left( x, y + \Delta y \right) + U \left( x, y - \Delta y \right) - 2 U \left( x, y \right)}{\Delta y^2} = -4 \pi \rho \left( x, y \right) \end{array} \end{equation}

If the grid spacing is the same in both the \( x \) and \( y \) coordinates, that is \( \Delta x = \Delta y \), then the equation reduces to

\begin{equation} \begin{array}{l} U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) + U \left( x, y + \Delta \right) \\ \hspace{2 in} + \hspace{1 mm} U \left( x, y - \Delta \right) - 4 U \left( x, y \right) = -4 \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Algebraically, the solution \( U \left( x, y \right) \) can be shown as

\begin{equation} \begin{array}{l} U \left( x, y \right) \simeq \dfrac{1}{4} \left[ U \left( x + \Delta, y \right) + U \left( x - \Delta, y \right) \right. \\ \hspace{2 in} \left. + \hspace{1 mm} U \left( x, y + \Delta \right) + U \left( x, y - \Delta \right) \right] + \pi \rho \left( x, y \right) \Delta^2 \end{array} \end{equation}

Lattice spacings can be identified as \( x = x_0 \pm \Delta \) and \( y = y_0 \pm \Delta \) which correspond to \(i \pm 1 \) and \( j \pm 1 \) where \( i, j = 0, 1, \ldots, N_{max - 1} \).  Thus, the finite difference algorithm appears as

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) + \pi \rho_{i, j}\Delta^2 \end{equation}

However, for simplicity Landau et al just ask to solve for the homogeneous equation

\begin{equation} U_{i, j} = \dfrac{1}{4} \left( U_{i + 1, j} + U_{i - 1, j} + U_{i, j + 1} + U_{i, j - 1} \right) \end{equation}

A few notes here from Landau et al.  For one a direct solution cannot be obtain from this algorithm.  That is, a process must be repeated several times in order for a solution to converge!  Thus, this is where efficiency and proficiency becomes very important in numerical methods.

First, a guess is imposed then swept over the domain which finds values for all variables at each point or node.  The solution has "converged" when the current generation does not change compared to the previous generation of values within a "tolerance" or error/precision (absolute value of the difference between old values and new) specified by the user or technique.  Another outcome could be that the values emerge as not expected which at that point the user must go back and figure out from where the error is coming.  Landau et al give the definition of relaxed as when the first guess has converged into a solution.

Landau et al pose "Does it always converge, and if so, does it converge fast enough to be useful?"

Another important note is the coordinate system used.  This problem was done in rectangular Cartesian coordinates.  Other problems may work better in other coordinate systems such as cylindrical or spherical depending upon the geometry of the problem in which case the Laplacian expansion \( \Delta^2 \) will be different.

One way to solve the algorithm is to utilize the Jacobi method.  This most basic method simply evaluates each point, sweeping across the domain.  The initialization and BCs symmetry are preserved in this way.  However, in order to speed up this method, a simple modification can be implemented where values that are known from the current run compared to other nodes can be updated immediately.  Thus, reducing the amount of needed storage by half since two nodes will be know for one in this case and the previous sweep will not be necessary to store since information is updated instantaneously when available.  This technique is called the Gauss-Seidel method.  The GS method leads to quicker convergence, less round-off error, and less needed storage.  However, symmetry is not conserved which the user hopes does not affected the solution outcome.

If the sweep begins at the top left corner then the GS algorithm looks like

\begin{equation} U_{i, j} = \dfrac{1}{4} \left[ U^{(old)}_{i + 1, j} + U^{(new)}_{i - 1, j} + U^{(old)}_{i, j + 1} + U^{(new)}_{i, j - 1} \right] \end{equation}

Octave code:

%This code was written by Tim Barber for the numerical PDE class taught by Dr. Parigger in the Summer of 2011

%This file corresponds specifically to HW set #1

%Problem 1.3 Square Wire Problem: Chapter 17 of the Landau-course book

discusses the square-wire problem, viz. elliptic PDE. Section 17.4.2

elaborates on the solution using the finite difference method. Your task

here is to (i) derive the finite difference algorithm Eq. 17.27, (ii)

implement the code using source programing (see the listed source

code), and (iii) display the results using gnuplot (along the lines of

Section 17.4.2).

clf

clc

clear

%input data manually

N_max = 100;

iter = 1000;

%initialize matrix

for i = 1:N_max

for j = 1:N_max

if j = 100

U(i, 100) = 99;

else

U(i, j) = 0;

j = j + 1;

end

end

i = i +1;

end

%iterate matrix using the Gauss-Seidel technique

for n = 1:iter + 1

for i = 2:N_max - 1

for j = 2:N_max - 1

U(i, j) = (1/4)*(U(i + 1, j) + U(i - 1, j)

+ U(i, j + 1) + U(i, j - 1));

j = j + 1;

end

i = i +1;

end

n = n + 1;

end

%plot of U(x, y)

mesh (U)

%label title and axes

title('HW 1.3', "fontsize", 20)

ylabel('y', "fontsize", 16)

xlabel('x', "fontsize", 16)

zlabel('U(x, y)', "fontsize", 16)

%commands to save plot as a .png file

print("HW_1_3_mesh.png")

replot

---

References:

[1] R. H. Landau, M. J. Páez, and C. C. Bordeianu. A Survey of Computational Physics -

a Introductory Computational Science, Princeton University Press, Princeton, New Jersey.

2008

Numerical Math - Stiffness

In order to answer the question of stiff ODEs, let's look at an example in A First Course in the Numerical Analysis of Differential Equations by Iserles [1].

Hoffman [2] states that stiffness appears for single linear and nonlinear ODEs, higher-order linear and nonlinear ODEs, and systems of linear and nonlinear ODEs.  He lists several definitions for stiffness:

1. An ODE is stiff if the step size required for stability is much smaller than the step size required for accuracy.

2. An ODE is stiff if it contains some components of the solution that decay rapidly compared to other components of the solution.

3. A system of ODEs is stiff if at least one eigenvalue of the system is negative and large compared to the other eigen values of the system.

4. From a practical point of view, an ODE is stiff if the step sixe based on cost (i.e., computational time) is too large to obtain an accurate (i.e., stable) solution.

Hoffman [2] also gives a couple of examples.

In progress...to be continued.

---

References:

[1] A. Iserles. A First Course in the Numerical Analysis of Differential Equations. 2nd ed. Cambridge University Press, Cambridge, UK. 2009

[2] J. D. Hoffman. Numerical Methods for Engineers and Scientists. 2nd ed. Marcel Dekker, Inc., New York, NY. 2001

Numerical Math - Adams methods

According to Iserles [1],

In progress...to be continued.

---

[1] A. Iserles. A First Course in the Numerical Analysis of Differential Equations. (2nd ed.).

Cambridge University Press, Cambridge, UK. 2009

Numerical Mathematics - Solving the heat diffusion equation

This post is influenced by a homework problem from a numerical PDE class I took this past summer (2011). This problem examines the DMOLCH routine and Gear's method for a heat diffusion equation. We are to reproduce data found in a journal paper.  The numerics to obtain can be found in the journal article "Temperature distribution in dental tissue after interaction with femtosecond laser pulses" by Pike et al in Applied Optics, Vol. 46,  Num. 34, Dec. 2007.  I also include some discussion on IMSL, FORTRAN, MOL and many other topics associated with this problem.  I will also have individual posts on these topics as well.

---

IMSL stands for the International Mathematics and Statistics Library and consist of numerous commercial software library utilized in numerical analysis in various programming languages organized by Rogue Wave Software [1, 2]. According to the 2003 Visual Numerics IMSL Fortran Library User’s Guide Math Library Volume 2 of 2 - Mathematical Functions in Fortran [3], the package routine DMOLCH is the double precision version of the routine MOLCH in the Fortran 77 interface which of course is single precision. In the Fortran 90 interface the specific interface names are S MOLCH and D MOLCH for double and single precision, respectively. The routine MOLCH according to the 2003 guide,

“Solves a system of partial differential equations of the form \( u_t = f \left( x, t, u, u_x, u_{xx} \right) \)  using the method of lines. The solution is represented with cubic Hermite polynomials.”

According to the 2010 guide, Visual Numerics IMSL Fortran Numerical Library User’s Guide Math Library Version 7.0 [4], the MOLCH routine has been “deprecated” by the new routine MMOLCH, where the definition remains the same as the quote above and the specific F90 interface names for single and double precision are S MMOLCH and D MMOLCH, respectivley.

The method of lines (MOL) is described by Wouwer et al [5] in the second section of the first chapter in the book entitled Adaptive Method of Lines [6]. The MOL attacks and aims to solve partial differential equations (PDEs) with the form \begin{equation} \textbf{u}_t = \textbf{f}(\textbf{u}), \quad \quad \textbf{x}_l < \textbf{x} < \textbf{x}_r, \quad \quad t > 0 \end{equation} where \begin{align*} \textbf{u}_t &= \dfrac{\partial \textbf{u}}{\partial t} \\ \textbf{u} &= \text{vector of dependent variables} \\ t &= \text{initial value independent variable} \\ \textbf{x}& = \text{boundary value independent variables} \\ \textbf{f} &= \text{spatial differential operator} = \textbf{f} \left( \textbf{x}, t, \textbf{u}, \textbf{u}_x,\textbf{u}_{xx}, \ldots \right) \end{align*} Wouwer et al [5] divide the MOL process into two sections. The first step discretizes the spatial derivatives, \( \textbf{u}_x, \textbf{u}_xx, \ldots \) using finite difference (FD) or finite element (FE) methods. The second step involves integrating the resulting system of semi-discrete ordinary differential equations (ODEs) in the initial value variable wrt to time, \( t \).

As an example of the MOL, Wouwer et al [5] utilize the linear advection equation to demonstrate the process. \[ \dfrac{\partial u}{\partial t} = -v \dfrac{\partial u}{\partial x} \] To begin, the equation is discretized in the spatial realm of \( x \) into an \( N \) number of grid points which are spaced evenly by the distance \( \delta x \). The spatial derivative, \( \partial u/ \partial x \), is then approximated by its second order centered FD beginning at grid point \( i \) up to the maximum discretization point \( N \) so that \[ \dfrac{\partial u}{\partial x} = -v \dfrac{u_{i+1} - u_{i-1}}{2 \Delta x} + O \left( \Delta x^2 \right), \quad i = 1, 2, \cdots, N \] Then assuming a value of \( v = 1 \) and substituting into \( \partial u/\partial t = -v \partial u/\partial x \) results in the equation \[ \dfrac{\mathrm{d} u}{\mathrm{d} t} = -\dfrac{u_{i+1} - u_{i-1}}{2 \Delta x}, \quad i = 1, 2, \cdots, N \]

Wouwer et al [5] next mention that at this point a system of \( N \) ODEs are generated which can then be integrated by an already developed by a library ODE integrator. An initial condition must then be specified next at the \( N \) grid points and the false points \( u_0 \) and \( u_{N + 1} \) for \( i = 1 \) and \( i = N \).  Thus, if an example of the equation is solved for two step functions separated by a time interval of \( t = 50 \)  or simply a square pulse, with a spatial length of \( N =101 \) or 100 intervals with a length of \( \Delta x \), the solution exhibits oscillatory behavior as seen in the figure by Wouwer et al [5]. This oscillation is of course false and is caused by the numerics of the problem known as a numerical distortion.

Wouwer et al [5] further discuss some issues here about this problem that I feel like are good points to intake. First is a continuation of the numerical oscillation problem. Wouwer et al express the problem that the oscillations do not dissipate with an increase in grid resolution in the spatial direction \( x \) or, in other words, and increase in the number of intervals, even though the FD equation is of the second order so you would expect as \( \Delta x \)  decreased the error would decrease \( O \left( \Delta x^2 \right) \).  Thus, the observation is made that will a reasonable method of approximation of a linear PDE broken down into a system of ODEs is accomplished the solutions are inaccurate.

In fact, Wouwer et al [5] explain that the accuracy of the numerical solution is dependent upon the smoothness of the analytical solution. For instance, rather than the two discontinuous jumps encountered with the square wave, a “more smooth” function would be a triangular pulse where \( u\left(x,t\right) \) is continuous but \( \partial u/\partial x \) is not continuous. As shown in the figure the MOL numerical solution (dots) matches much closer with the actual solution (solid line). Finally, if  the initial condition is even “smoother” such as a cosine function, then the oscillations diminish even more.  Thus, if the initial condition contains non-smooth characteristics such as discontinuities, then the discontinuities spread through the time and space variables which is indicative of hyperbolic PDEs.

However, Gear’s method is also utilized by Pike et al [8]. Gear’s method as defined by Hoffman [7] was developed by Gear in 1971 [9]. Gear’s method contains much larger stability limits for a series of implicit FD equations (FDEs). The formula appears as

\[ y_{n+1} = \gamma \left[ \beta h f_{n+1} + \left( \alpha_0 y_n + \alpha_{-1} y_{n-1} + \alpha_{-2} y_{n-2} + \cdots \right)\right] \]

where k is the global order of the FDE and the coefficients \( \gamma \), \( \beta \), \( \alpha_i \left( i = 0, 1, 2, \ldots \right) \) are shown in the table  from Hoffman [7].

In addition, an interesting comment from Hoffman [7] cites that the Gear FDEs are in a FORTRAN package named LSODE which came about from researchers at Lawrence Livermore National Laboratory. According to Hoffman, this package utilizes

"...an elaborate error control procedure based on using various-order Gear FDEs in conjunction with step size halving and doubling."

Additional notes from Hoffman explain that stiff ODEs contain extra issues for solutions.  Usually explicit methods do not work for stiff ODEs.  However, implicit methods such as the Gear method do work for stiff problems.  Hoffman defines stiff ODEs as where

"...at least one eigenvalue has a large negative real part which causes the corresponding component of the solution to vary rapidly compared to the typical scale of variation displayed by the rest of the solution."

Velten [10] states that stiff ODEs pose issues since the numerical solutions may oscillate as demonstrated beforehand. Thus, small step sizes must be implemented to obtain a correct solution. Also of note is Velten’s mention of a few software packages such as Maxima and R.  Velten gives another name for the Gear method as the backward differentiation formulas (BDF) method. A part of R’s odesolve package is lsoda which stands for Livermore Solver for Ordinary Differential Equations (automatic) which came about from Petzold and Hindmarsh at California’s Lawrence Livermore National Laboratory [11, 12]. The lsoda package has a unique feature in which it automatically switches between methods for stiff and nonstiff sets of ODEs which is one of the main differences from lsode. As a further side note, according to Soetaert et al [13] the package odesolve is planned to be deprecated and deSolve is the package successor. The nonstiff case utilizes the Adams method with a variable stepsize and a variable order of up to 12th order. For the stiff case of ODEs, the lsoda package utilizes the BDF method with a variable stepsize and a variable order of up to 5th order.

In order to answer the question of stiff ODEs, let's look at an example in A First Course in the Numerical Analysis of Differential Equations by Iserles [14].

Hoffman [7] states that stiffness appears for single linear and nonlinear ODEs, higher-order linear and nonlinear ODEs, and systems of linear and nonlinear ODEs.  He lists several definitions for stiffness:

1. An ODE is stiff if the step size required for stability is much smaller than the step size required for accuracy.

2. An ODE is stiff if it contains some components of the solution that decay rapidly compared to other components of the solution.

3. A system of ODEs is stiff if at least one eigenvalue of the system is negative and large compared to the other eigen values of the system.

4. From a practical point of view, an ODE is stiff if the step sixe based on cost (i.e., computational time) is too large to obtain an accurate (i.e., stable) solution.

Hoffman [7] also gives a couple of examples.

Also, since it is mentioned let’s talk about Adams methods. According to Iserles [14],


In progress...to be continued.

---

References:

[1] IMSL Numerical Libraries. 07-23-2011, http://en.wikipedia.org/wiki/IMSL Numerical Libraries

[2] IMSL R Numerical Libraries. imsl-numerical-libraries.aspx 07-23-2011, http://www.roguewave.com/products/

[3] Visual Numerics IMSL Fortran Library User’s Guide Math Library Volume 2 of 2 - Mathematical Functions in Fortran. 07-23-2011, http://www.absoft.com/Support/Documentation/MathV2.pdf, 2003

[4] Visual Numerics IMSL Fortran Numerical Library User’s Guide Math Library Version 7.0. 07-23-2011, http://www.roguewave.com/portals/0/products/imsl-numerical-libraries/fortran-library/docs/7.0/math/math.pdf, 2010

[5] A. Vande Wouwer, P. Saucez, & W. E. Schiesser. Chapter 1 Introduction. A. Vande Wouwer, P. Saucez, & W. E. Schiesser (editors). Adaptive Method of Lines. Chapman & Hall/CRC, Boca Raton, FL. 2001

[6] A. Vande Wouwer, P. Saucez, & W. E. Schiesser (editors). Adaptive Method of Lines. Chapman & Hall/CRC, Boca Raton, FL. 2001

[7] J. D. Hoffman. Numerical Methods for Engineers and Scientists. 2 edition. Marcel Dekker, Inc., New York, NY. 2001

[8] P. Pike, C. Parigger, R. Splinter, and P. Lockhart. Temperature distribution in dental tissue after interaction with femtosecond laser pulses. Applied Optics, Volume 46, Number 34, December 2007

[9] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, New Jersey. 1971

[10] K. Velten. Mathematical Modeling and Simulation: Introduction for Scientists and Engineers, Wiley-VCH, Weinheim, Germany. 2009

[11] A. C. Hindmarsh. ODEPACK, a systematized collection of ODE solvers, Scientific Computing, North-Holland, Amsterdam, The Netherlands, pp. 55–64. 1983

[12] L. R. Petzold. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM Journal on Scientific and Statistical Computing, 4, 136–48. 1983

[13] K. Soetaert, T. Petzoldt, and R. W. Setzer. Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33, 4, February 2010


[14] A. Iserles. A First Course in the Numerical Analysis of Differential Equations. (2nd ed.).
Cambridge University Press, Cambridge, UK. 2009

Additional keywords:

finite difference equations (FDEs), Gear's method, LSODE, LSODA, Adams method, deSolve, R, Maxima, backward differentiation formulas (BDF) method, International Mathematics and Statistics Library (IMSL), Rogue Wave Software, DMOLCH, MOLCH, Fortran, Fortran 90, MMOLCH, Fortran 77, S MOLCH, D MOLCH, S MMOLCH, D MMOLCH, Method of Lines (MOL), Partial Differential Equations (PDEs), Finite Difference (FD), Finite Element (FE), ordinary differential equations (ODEs)