473,224 Members | 1,368 Online

# numpy help

Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray.linear_algebra.inverse and
numpy.oldnumeric.linear_algebra.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).
This is the error message I receive:

Traceback (most recent call last):
File "C:\Documents and Settings\Chris &
Esther\Desktop\636_hw5_2\elen636_hw5_2.py", line 60, in <module>
matrix_inverse =
numpy.numarray.linear_algebra.generalized_inverse( matrix)
File
"C:\Python25\lib\site-packages\numpy\oldnumeric\linear_algebra.py", line
59, in generalized_inverse
return linalg.pinv(a, rcond)
File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line
557, in pinv
u, s, vt = svd(a, 0)
File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line
485, in svd
a = _fastCopyAndTranspose(t, a)
File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line
107, in _fastCopyAndTranspose
cast_arrays = cast_arrays + (_fastCT(a.astype(type)),)
TypeError: can't convert complex to float; use abs(z)

I've tried inverting small complex matrices and it worked fine. Does
anyone know why it won't work for this larger matrix? Any ideas how I
can work around this problem and get the correct inverse matrix?

Chris

P.S. elen636_math.py is my personal library of functions I've create to
solve the problem while elen636_hw5_2.py is the program that I'm
actually running
# Purpose:
# This is a library of functions for ELEN 636 that
# so far has the ability to calculate the Sine and
# Cosine integrals as well as the mutual impedance
# between two parallel antennas.
#
# Author: Christopher Smith
# E-mail: mc*****@tamu.edu
# Date: 10/30/2006

###############
### NOTE: The functions below for the sine and cosine integrals are similar
### to the functions I turned in for homework assignment 4 problem 6
### except that I added the ability to check for convergence.
### I also added the factor into the mutual impedance formula so that the
### answer is given in terms of the terminal input impedance instead of
### the loop impedance as it was formerly giving.
###############

# depends on the math library
from math import *
import numpy.numarray, numpy

def factorial(n):
"""
This function calculates the factorial of a number.
"""
sum = 1.0
for m in range(1, int(n)+1):
sum = float(m)*sum
return sum

def Si(x):
"""
This function computes the sine integral. It uses a power series
expansion that can be found in Abramowitz and Stegun's math
functions reference book.
"""
start = 0.0
stop = 10.0
sine_int = 0.0
convergence = 1.0*10**(-6) # want to have the difference between
# the last run and this run below
# this value
while 1:
for n in range(int(start), int(stop)):
n = float(n)
sine_int += ((-1)**n)*x**(2*n +1)/((2*n+1)*factorial(2*n+1))
sine_int_new = sine_int + ((-1.)**stop)*x**(2.*stop +1.)/((2.*stop+1.)*factorial(2.*stop+1.))
converge_check = sine_int_new - sine_int
if abs(converge_check) < convergence:
break
else:
start = stop
stop += 5.0
return sine_int_new

def Ci(x):
"""
This function computes the cosine integral. It uses a power series
expansion that can be found in Abramowitz and Stegun's math
functions reference book.
"""
start = 1.0
stop = 10.0
convergence = 1.0*10.**(-6) # want to have the difference between
# the last run and this run below
# this value
# The first number in the sum is Euler's constant to 10 digits
cosine_int = 0.5772156649 + log(x)
while 1:
for n in range(int(start), int(stop)):
m = float(n)
cosine_int = cosine_int +((-1)**m)*x**(2*m)/((2*m)*factorial(2*m))
cosine_int_new = cosine_int + ((-1)**stop)*x**(2*stop)/((2*stop)*factorial(2*stop))
converge_check = cosine_int_new - cosine_int
if abs(converge_check) < convergence:
break
else:
start = stop
stop += 5.0
#print stop
return cosine_int_new
def mutual_impedance(length1_tot, length2_tot, stagger, d):
"""
This function computes the mutual impedance between two antennas
for the Parallel in Echelon Configuration. The formulas are taken
from a paper by Howard King, "Mutual Impedance of Unequal Length
Antennas in Echelon"

NOTE: all measurements should be entered in wavelengths
"""

# stagger (this is the vertical separation between antenna centers)
# d (this is the horizontal separation between the antennas)

# length1 and length2 are the half length of the antennas, this is
# to conform to King's formulas
length1 = length1_tot/2.0
length2 = length2_tot/2.0

# vertical separation between center of antenna 1 and bottom of antenna 2
h = stagger - length2

# wave propagation constant
beta = 2.0*pi

# formulas to put into mutual impedance equation
u0 = beta*(sqrt(d**2 +(h -length1)**2) +(h -length1))
v0 = beta*(sqrt(d**2 +(h -length1)**2) -(h -length1))
u0prime = beta*(sqrt(d**2 +(h +length1)**2) -(h +length1))
v0prime = beta*(sqrt(d**2 +(h +length1)**2) +(h +length1))
u1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) +(h -length1 +length2))
v1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) -(h -length1 +length2))
u2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) -(h +length1 +length2))
v2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) +(h +length1 +length2))
u3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)**2) +(h -length1 +2.0*length2))
v3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)**2) -(h -length1 +2.0*length2))
u4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)**2) -(h +length1 +2.0*length2))
v4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)**2) +(h +length1 +2.0*length2))
w1 = beta*(sqrt(d**2 +h**2) -h)
y1 = beta*(sqrt(d**2 +h**2) +h)
w2 = beta*(sqrt(d**2 +(h +length2)**2) -(h +length2))
y2 = beta*(sqrt(d**2 +(h +length2)**2) +(h +length2))
w3 = beta*(sqrt(d**2 +(h +2.0*length2)**2) -(h +2.0*length2))
y3 = beta*(sqrt(d**2 +(h +2.0*length2)**2) +(h +2.0*length2))
#print u0,v0,u0prime,v0prime,u1,v1,u2,v2,u3,v3,u4,v4,w1,y 1,w2,y2,w3,y3

# real part of the mutual impedance between two antennas
R12 = 15*(cos(beta*(length1 -h))*(Ci(u0) +Ci(v0) -Ci(u1) -Ci(v1)) \
+sin(beta*(length1 -h))*(-Si(u0) +Si(v0) +Si(u1) -Si(v1)) \
+cos(beta*(length1 +h))*(Ci(u0prime) +Ci(v0prime) -Ci(u2) -Ci(v2)) \
+sin(beta*(length1 +h))*(-Si(u0prime) +Si(v0prime) +Si(u2) -Si(v2)) \
+cos(beta*(length1 -2.0*length2 -h))*(-Ci(u1) -Ci(v1) +Ci(u3) +Ci(v3)) \
+sin(beta*(length1 -2.0*length2 -h))*(Si(u1) -Si(v1) -Si(u3) +Si(v3)) \
+cos(beta*(length1 +2.0*length2 +h))*(-Ci(u2) -Ci(v2) +Ci(u4) +Ci(v4)) \
+sin(beta*(length1 +2.0*length2 +h))*(Si(u2) -Si(v2) -Si(u4) +Si(v4)) \
+2.0*cos(beta*length1)*cos(beta*h)*(-Ci(w1) -Ci(y1) +Ci(w2) +Ci(y2)) \
+2.0*cos(beta*length1)*sin(beta*h)*(Si(w1) -Si(y1) -Si(w2) +Si(y2)) \
+2.0*cos(beta*length1)*cos(beta*(2.0*length2 +h))*(Ci(w2) +Ci(y2) -Ci(w3) -Ci(y3)) \
+2.0*cos(beta*length1)*sin(beta*(2.0*length2 +h))*(-Si(w2) +Si(y2) +Si(w3) -Si(y3)))

# imaginary part of the mutual impedance between two antennas
X12 = 15*(cos(beta*(length1 -h))*(-Si(u0) -Si(v0) +Si(u1) +Si(v1)) \
+sin(beta*(length1 -h))*(-Ci(u0) +Ci(v0) +Ci(u1) -Ci(v1)) \
+cos(beta*(length1 +h))*(-Si(u0prime) -Si(v0prime) +Si(u2) +Si(v2)) \
+sin(beta*(length1 +h))*(-Ci(u0prime) +Ci(v0prime) +Ci(u2) -Ci(v2)) \
+cos(beta*(length1 -2.0*length2 -h))*(Si(u1) +Si(v1) -Si(u3) -Si(v3)) \
+sin(beta*(length1 -2.0*length2 -h))*(Ci(u1) -Ci(v1) -Ci(u3) +Ci(v3)) \
+cos(beta*(length1 +2.0*length2 +h))*(Si(u2) +Si(v2) -Si(u4) -Si(v4)) \
+sin(beta*(length1 +2.0*length2 +h))*(Ci(u2) -Ci(v2) -Ci(u4) +Ci(v4)) \
+2.0*cos(beta*length1)*cos(beta*h)*(Si(w1) +Si(y1) -Si(w2) -Si(y2)) \
+2.0*cos(beta*length1)*sin(beta*h)*(Ci(w1) -Ci(y1) -Ci(w2) +Ci(y2)) \
+2.0*cos(beta*length1)*cos(beta*(2.0*length2 +h))*(-Si(w2) -Si(y2) +Si(w3) +Si(y3)) \
+2.0*cos(beta*length1)*sin(beta*(2.0*length2 +h))*(-Ci(w2) +Ci(y2) +Ci(w3) -Ci(y3)))

R12_in = R12/(sin(beta*length1)*sin(beta*length2))
X12_in = X12/(sin(beta*length1)*sin(beta*length2))
impedance = (R12_in, X12_in)
return impedance
def top_row_matrix(length1, length2, stagger, stagger_image, radius, m, n):
"""
This function will find the top row of a mutual impdedance matrix
over a ground plane. From the top row we can find the overall matrix
since it is a block Toeplitz matrix
"""

z = [] # list to store our impedance values

# index to step over for the staggering and separation between antennas
stagger_range = range(0, m)
separation_range = range(0, n)

# calculate the mutual impedance values for the real planar array
# the first loop gives us the stagger between rows while the second
# loop gives us the separation between dipoles on the same row
for m in stagger_range:
for n in separation_range:
h = stagger*m
d = stagger*n
if d == 0: d = radius
trans = mutual_impedance(length1, length2, h, d)
z.append(complex(trans[0],trans[1]))

# Suppose the real antenna array is a plane in x-y at z = 0 and the
# imaginary antenna array is a plane in x-y at z = some spacing
# Since the mutual impedance would include spacing in the
# x,y and z directions we need to condense it down to two spacings
# so that we can input the spacings into our formulas. For the
# spacing between dipoles we will use the hypotenuse formed by the
# x and z coordinates. The y coordinate will be the staggering between
# dipoles.
for m in stagger_range:
for n in separation_range:
separation = sqrt((stagger_image)**2+(stagger*n)**2)
trans = mutual_impedance(length1, length2, stagger*m, separation)
z.append(complex(trans[0],trans[1]))

return z
def Toeplitz(a):
"""
This function takes a list in. The list represents
the top row in a Toeplitz matrix. From this row we
can fill the rest of the matrix and then output it.

Note: This function has a dependency on the numpy
library.
"""

# need the length of the list being input
n = len(a)

# fill a matrix the size we need with zeros
matrix = numpy.zeros((n,n), dtype = type(a))

a_new = []
for i in range(-n+1,n):
a_new.append( a[ abs(i) ] )

for i in range(0,n):
matrix[i,:] = a_new[ (n-1) -i : (2*n -1) -i ]

return matrix
def BlockToeplitz(a, M, N):
"""
This function takes a list of N Toeplitz matrices in
and creates a Block Toeplitz matrix from them.

Note: This function has a dependency on the numpy
library.
This function is also specifically for an array
over a groundplane problem. Which means that
the matrix is twice the length and width of a
normal array matrix. The code would need to be
modified to apply to a normal Block-Toeplitz.
"""
# fill a matrix the size we need with zeros
matrix = numpy.zeros((2*M*N,2*M*N), dtype = type(a))

a_new = []
for i in range(-2*N+1,2*N):
a_new.append( a[ abs(i) ] )

for i in range(0, 2*N):
i_new = i*M
for j in range(0, 2*N):
j_new = j*M
matrix[ i_new:i_new+N, j_new:j_new+N ] = a_new[2*N-1-j+i]

return matrix
def h_array_over_gndplane(M, N):
"""
This function generates the voltage excitation matrix for an
array over a ground plane. The antennas are parallel to the
ground plane which means that according to image theory the
current along the antenna images are in the opposite direction
from those of the real array.

Note: This function has a dependency on the numpy library.
"""

h = []
for n in range(0, 2*M*N):
if n < M*N:
h.append(1.0)
else:
h.append(-1.0)
h = numpy.numarray.array((h), shape = (2*M*N,1))
return h
# Purpose:
# This program will calculate the impedance matrix
# for a planar antenna array over a ground plane.
# This is for ELEN 636 homework 5 problem 2
#
# Author: Christopher Smith
# E-mail: mc*****@tamu.edu
# Date: 10/30/06

# library so arrays (matrices) can be used
import numpy.numarray, numpy
# library so we can perform linear algebra operations
# like inverse( )
import numpy.numarray.linear_algebra
# import special functions
from elen636_math import mutual_impedance, top_row_matrix, Toeplitz
from elen636_math import BlockToeplitz, h_array_over_gndplane
# math library
from math import *
# length and radius of dipoles (wvlgths)
# horizontal (same in x and y) stagger of matrix (wvlgths)
# vertical separation between array and it's image (wvlgths)
# or twice the vertical separation between array and ground plane
# and the number of elements in the y and x directions m, n
length1 = 0.5
length2 = length1
radius = 0.0025
stagger = 0.55
stag_image = 0.5
M = 7
N = 7

# function to give us back the top row of the m x n matrix
z = top_row_matrix(length1, length2, stagger, stag_image, radius, M, N)

# We now have a list with 2*M*N elements. Using symmetry we can create
# N NxN arrays which correspond to N Toeplitz matrices for our N x N
# Block Toeplitz matrix of 2*M*N x 2*M*N elements. So from those N
# matrices we can completely fill the matrix.
toeplitz_matrices = [] # a list of N Toeplitz matrices
for i in range(0, 2*M*N, N):
toeplitz_matrices.append( Toeplitz( z[i:i +N] ) )

# Now we fill our Block-Toeplitz matrix
matrix = BlockToeplitz(toeplitz_matrices, M, N)

# generate the voltage excitations of all the dipoles
h = h_array_over_gndplane(M, N)

# now we compute Cn by multiply the inverse mutual impedance matrix
# and multiplying it by the h matrix
matrix_inverse = numpy.numarray.linear_algebra.inverse(matrix)
Cn = numpy.numarray.matrixmultiply(matrix_inverse,h)

input_impedance_matrix = h/Cn

#input_impedance_element25 = input_impedance_matrix[24]

# creates and opens files to output results
try:
file = open( 'mutual_impedance_results.txt', "w" )
except IOError, message:
print >sys.stderr, "File could not be opened:" , message
sys.exit(1)

print >file, input_impedance_matrix
file.close()

_______________________________________________
Tutor maillist - Tu***@python.org
http://mail.python.org/mailman/listinfo/tutor
Nov 3 '06 #1
2 3917
Chris Smith wrote:
Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray.linear_algebra.inverse and
numpy.oldnumeric.linear_algebra.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).
try scipy.linalg.inv etc.

(you are using outdated old NumPy modules with recent array type which sometimes creates type conflicts.

robert

Nov 3 '06 #2
robert wrote:
Chris Smith wrote:
>>Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray.linear_algebra.inverse and
numpy.oldnumeric.linear_algebra.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).

try scipy.linalg.inv etc.

(you are using outdated old NumPy modules with recent array type which sometimes creates type conflicts.

robert
I tried that one also, and it didn't work.

Chris

Nov 3 '06 #3

This thread has been closed and replies have been disabled. Please start a new discussion.

### Similar topics

 2 by: Justin Lemkul | last post by: Hello all, I am hoping someone out there will be able to help me. I am trying to install a program that utilizes NumPy. In installing NumPy, I realized that I was lacking Atlas. I ran into... 1 by: drife | last post by: Hello, I use the Python Numeric package extensively, and had been an avid user of the "old" scipy. In my view, both pieces of software are truly first rate, and have greatly improved my... 13 by: Gary Wessle | last post by: Hi I am trying to install NumPy in my debian/testing linux 2.6.15-1-686. with no numpy for debian/testing, I am left alone, since the experimental version available by debian will result in a... 15 by: greg.landrum | last post by: After using numeric for almost ten years, I decided to attempt to switch a large codebase (python and C++) to using numpy. Here's are some comments about how that went. - The code to... 2 by: robert | last post by: in Gnuplot (Gnuplot.utils) the input array will be converted to a Numeric float array as shown below. When I insert a numpy array into Gnuplot like that below, numbers 7.44 are cast to 7.0 Why is... 1 by: Jianzhong Liu | last post by: Hello, Guys, I have a question about the linear_least_squares in Numpy. My linear_least_squares cannot give me the results. I use Numpy1.0. The newest version. So I checked online and get... 3 by: Duncan Smith | last post by: Hello, Since moving to numpy I've had a few problems with my existing code. It basically revolves around the numpy scalar types. e.g. ------------------------------------------------ array(,... 5 by: Marc Oldenhof | last post by: Hello all, I'm pretty new to Python, but use it a lot lately. I'm getting a crazy error trying to do operations on a string list after importing numpy. Minimal example: Python 2.5.1... 3 by: Sean Davis | last post by: I have a set of numpy arrays which I would like to save to a gzip file. Here is an example without gzip: b=numpy.ones(1000000,dtype=numpy.uint8) a=numpy.zeros(1000000,dtype=numpy.uint8) fd =... 1 by: isladogs | last post by: The next online meeting of the Access Europe User Group will be on Wednesday 6 Dec 2023 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, Mike... 0 by: veera ravala | last post by: ServiceNow is a powerful cloud-based platform that offers a wide range of services to help organizations manage their workflows, operations, and IT services more efficiently. At its core, ServiceNow... 3 by: isladogs | last post by: The next Access Europe meeting will be on Wednesday 3 Jan 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). For other local times, please check World Time Buddy In... 0 by: jianzs | last post by: Introduction Cloud-native applications are conventionally identified as those designed and nurtured on cloud infrastructure. Such applications, rooted in cloud technologies, skillfully benefit from... 0 by: mar23 | last post by: Here's the situation. I have a form called frmDiceInventory with subform called subfrmDice. The subform's control source is linked to a query called qryDiceInventory. I've been trying to pick up the... 2 by: jimatqsi | last post by: The boss wants the word "CONFIDENTIAL" overlaying certain reports. He wants it large, slanted across the page, on every page, very light gray, outlined letters, not block letters. I thought Word Art... 2 by: isladogs | last post by: The next Access Europe meeting will be on Wednesday 7 Feb 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:30 (7.30PM). In this month's session, the creator of the excellent VBE... 0 by: egorbl4 | last post by: Скачал я git, хотел начать настройку, а там вылезло вот это Что это? Что мне с этим делать? ... 0 by: MeoLessi9 | last post by: I have VirtualBox installed on Windows 11 and now I would like to install Kali on a virtual machine. However, on the official website, I see two options: "Installer images" and "Virtual machines"....

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.