# [Numpy-discussion] First Numerical Python code...comments?

Tyler Hayes hayes.tyler@gmail....
Sat Mar 3 17:05:01 CST 2007

```Hello All:

Well, I recently made the switch to Python for scientific computing
(using NumPy and Matplotlib/PyLab) and got my first program to work.
It is a simple 1D finite element code for a tightly stretched wire and
is a textbook example.

To solve it, I also created a symmetric, banded matrix solver module
from a Fortran90 routine using f2py (thank you Pearu for your help the
other day to get f2py working).

Anyways, my background is in Fortran and I was got the program to run,
but I can't help but wonder if I did not use the most efficient method
for calcualting the matrix values, etc.

I would appreciate your comments on my code that I have attached
below. My concern is that I used some "for" loops as I would have
using Fortran, but perhaps there is a more intelligent (i.e.,
efficient) way using native NumPy matirx routines?

I would love to hear what you think.

Also, while I'm at it. I have a simple question about importing NumPy,
SciPy, and PyLab. If I import PyLab, it says I also import NumPy. Do I
ever need to call in SciPy as well? Or is this the same thing?

Cheers,

t.

The "code"

#!/usr/bin/env python
#
# This is the example code from Thompson for a tight wire problem
#
# Coded by Tyler Hayes March 1, 2007
#
#################################################################

# Load the pylab module for plotting (it also imports NumPy)
from pylab import *
# Load the sparse matrix solver symgauss created with f2py
import symgauss as G
"""
This is the example program from Thompson which will plot the

DEPENDS:
This code requires the symgauss module to solve the sparse
symmetric matrix.
"""

# Initialize data
#----------------------------------------------------------------

# Parameters
tens  = 20000.0
numel = 8
numnp = numel + 1
IB    = 2
BLAST = 1.0e6

# Boundary conditions
npbc     = zeros(numnp,float)
npbc[0]  = 1.0
npbc[-1] = 1.0

# Displacement
U = [0.0]*numnp

# Element X-coordinates
xord     = zeros(numnp,float)
xord[1]  = 2.0
xord[2]  = 4.0
xord[3]  = 5.0
xord[4]  = 6.0
xord[5]  = 7.0
xord[6]  = 8.0
xord[7]  = 10.0
xord[8]  = 12.0

# Force data vector
F    = -1.0*ones(numnp,float)
F[1] = -2.0
F[2] = -1.5
F[-2]= -2.0
F[-3]= -1.5

# Initialize SK --- the sparse symmetrix stiffness matrix "K"
SK = zeros((numnp,IB),float)

# Main routine
#----------------------------------------------------------------

# Fill in SK with proper values
for i in range(numel):
RL = xord[i+1] - xord[i]
RK = tens/RL
SK[i][0] = SK[i][0] + RK
SK[i][1] = SK[i][1] - RK
IP1 = i+1
SK[IP1][0] = SK[IP1][0] + RK

# Set boundary conditions
for i in range(numnp):
if npbc[i] == 1:
SK[i][0] = SK[i][0]*BLAST
F[i] = U[i]*SK[i][0]

# Call Symmetric Sparse matrix solver
U = G.symgauss(SK,F,numnp,IB)

# Visualize the data
#----------------------------------------------------------------

# Plot data using matplotlib
plot(xord,U)
title('Deflection of a tighly stretched wire')
xlabel('Distance along wire (normalized)')
ylabel('Deflection of wire (normalized)')
grid(True)
show()

--
Tyler Joseph Hayes
600 Talbot St. -- Apt. 812
London, Ontario
N6A 5L9

Tel   : 519.435.0967
Fax  : 519.661.3198
Cell  : 416.655.7897
email: thayes@uwo.ca

GPG Key ID# 0x3AE05130
```