Chapter 01

Introduction

This Julia notebook accompanies the book "Numerical Methods for Engineers: by D. V. Griffiths and I.M. Smith (NMfE).

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

NMfE.jl is useful by itself. It can also be viewed as an introduction to the methods used in a follow on package I'm working on, CSoM.jl, which makes the Finite Element toolkit explained in "Programming the Finite Element Method" (PtFEM), also by D. V. Griffiths and I.M. Smith, available in Julia.

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 p_xx_yy_for.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 the ..._for.jl program, 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 p_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 such as fromSkyline().

The nmex example (briefly discussed in section 1.5, explained in chapter 4)



In [26]:

    
using NMfE, IterativeSolvers

For 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 [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.