Chapter 01

Introduction

These Julia notebooks accompany the book "Numerical Methods for Engineers" by D. V. Griffiths and I.M. Smith NMfE.

Neither the Julia package NMfE.jl nor these notebooks are a replacement for NMfE. The intention is to use the package or the notebooks as a replacement for the Fortran programs while studying the book.

NMfE by itself provides a practical, step by step introduction to the programming of a variety of numerical methods.

NMfE can also be viewed as an introduction to the methods used in a follow on book "Programming the Finite Element Method" PtFEM, also by D. V. Griffiths and I.M. Smith.

My intention is to develop a second package, PtFEM.jl, and associated notebooks which make the Finite Element toolkit explained in PtFEM available in Julia.

I have not started on that package in earnest, although I have experimented in a temporary package called CSoM.jl CSoM covering approximately the first six chapters of PtFEM. Once PtFEM.jl starts to take shape, I will deprecate CSoM.jl.

Prerequisite requirements

Using the notebooks requires a couple of preliminary steps:

Install Julia
Inside the REPL of Julia (do steps 1 to 6 just once):
1. Pkg.clone("https://github.com/goedman/NMfE.jl")
  
  and on e.g. OSX, to enable the notebooks:
2. ENV("PYTHON") = ""
3. ENV("JUPYTER") = ""
4. Pkg.add("IJulia") # or Pkg.build("Julia") is it is already installed
5. Pkg.add("PyPlot")
6. Pkg.add("Gadfly")
  
  Still inside the REPL:
7. using IJulia
8. @async notebook()
  
  This should start a browser and the Jupyter/Julia notebook interface.
  
  Some more info can be found at https://github.com/stevengj/julia-mit.

Approach

The approach taken in the NMfE.jl package is to provide 2 versions of most programs described in the NMfE book.

The first version is named f_xx_yy.jl and is pretty much a direct translation of the Fortran program (i.e. $\mathbf{Program}$ $\mathbf{xx.yy}$) in the book.

Often, at the end of f_xx_yy.jl programs, a Julia solution is shown and the results are compared using an @assert statement.

If deemed useful a 'pure' Julia version is also available, labeled j_xx_yy.jl.

Some of the Fortran functions that come with the book are contained in a separate library (nm_lib). Equivalent Julia functions are provided in the /src directory of the NMfE Julia package NMfE.jl.

Note

The Fortran programs that come with both books use the Skyline format to store sparse matrices. My intention is to use the Julia sparse matrices in all versions. I do provide a few conversion functions, e.g. fromSkyline().

The j_01_01 example (program nmex, a variant of p_04_01 from chapter 4)

Vector Iteration for "Largest" Eigenvalue and Eigenvector



In [26]:

    
using NMfE

In Julia, to solve an equation like:

$[\mathbf{A}]\{\mathbf{x}\}=\{\mathbf{b}\}$

we will provide $[\mathbf{A}]$ and intial values for $\{\mathbf{x_0}\}$:



In [27]:

    
A = [10. 5. 6.; 5. 20. 4.; 6. 4. 30.]









    Out[27]:





3x3 Array{Float64,2}:
 10.0   5.0   6.0
  5.0  20.0   4.0
  6.0   4.0  30.0

Initial guess for $\{\mathbf{x}\}$:



In [28]:

    
x0 = [1., 1., 1.];

A ';' at the end of the last line in a cell suppresses the output.



In [29]:

    
tol = 1.0e-5
limit = 100;

In this case I have opted for capturing the core part of program nmex in a Julia function nmex() stored in the /src/ch01 directory of the NMfE Julia package NMfE.jl:



In [30]:

    
(iters, eigval, eigvec) = nmex(A, x0; tol=tol, limit=limit, show=4);









    



1 3 40.0 [0.525,0.725,1.0]
2 3 36.05 [0.412621359223301,0.585991678224688,1.0]
3 3 34.81969486823856 [0.37496514638518225,0.5107149970125474,1.0]
4 3 34.292650866361285 [0.3587715192056842,0.46917124415013295,1.0]



In [31]:

    
println("No of iters:                   $iters")
println("Largest eigen value:           $eigval")
println("Corresponding eigen vector:    $eigvec")









    



No of iters:                   19
Largest eigen value:           33.70929144206253
Corresponding eigen vector:    [0.3001549575700916,0.3664458330546978,0.8806954370739895]

Julia has the built-in function eig():



In [32]:

    
?eig









    



search: eig eigs eigmin eigmax eigvec eigval eigvecs eigvals eigfact eigvals!







    Out[32]:




..  eig(A,[irange,][vl,][vu,][permute=true,][scale=true]) -> D, V

Computes eigenvalues and eigenvectors of ``A``. See :func:`eigfact` for
details on the ``balance`` keyword argument.

.. doctest::

   julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
   ([1.0,3.0,18.0],
   3x3 Array{Float64,2}:
    1.0  0.0  0.0
    0.0  1.0  0.0
    0.0  0.0  1.0)

``eig`` is a wrapper around :func:`eigfact`, extracting all parts of the
factorization to a tuple; where possible, using :func:`eigfact` is
recommended.

..  eig(A, B) -> D, V

Computes generalized eigenvalues and vectors of ``A`` with respect to ``B``.

``eig`` is a wrapper around :func:`eigfact`, extracting all parts of the
factorization to a tuple; where possible, using :func:`eigfact` is
recommended.



In [33]:

    
jeigval, jeigvec = eig(A)









    Out[33]:





([7.141760285928215,19.14906123147525,33.70917848259656],
3x3 Array{Float64,2}:
  0.933383  -0.196733  0.300153
 -0.303243  -0.87964   0.366438
 -0.191936   0.433046  0.880699)

No response to the '@assert ...' macro means 'true':



In [34]:

    
@assert round(maximum(jeigval), 3) == round(eigval, 3)



In [35]:

    
@assert round(maximum(jeigval), 4) == round(eigval, 4)









    



LoadError: AssertionError: round(maximum(jeigval),4) == round(eigval,4)
while loading In[35], in expression starting on line 1

The values agree up to the 4th decimal, hence above AssertionError.