Appending Numpy Array with Another Numpy Array

Question

I have a NumPy array with equations solved symbolically, with constants a and b. Here's an example of the cell at index (2,0) in my array "bounds_symbolic":

-a*sqrt(1/(a**6*b**2+1))

I also have an array, called "a_values", that I would like to substitute into my "bounds_symbolic" array. I also have the b-value set to 1, which I would also like to substitute in. Keeping the top row of the arrays intact would also be nice.

In other words, for the cell indexed at (2,0) in "bounds_symbolic", I want to substitute all of my a and b-values into the equation, while extending the column to contain the substituted equations. I then want to do this operation for the entirety of the "bounds_symbolic" array.

Here is the code that I have so far:

import sympy
import numpy as np

a, b, x, y = sympy.symbols("a b x y")

# Equation of the ellipse solved for y
ellipse = sympy.sqrt((b ** 2) * (1 - ((x ** 2) / (a ** 2))))

# Functions to be tested
test_functions = np.array(
    [(a * b * x), (((a * b) ** 2) * x), (((a * b) ** 3) * x), (((a * b) ** 4) * x), (((a * b) ** 5) * x)])

# Equating ellipse and test_functions so their intersection can be symbolically solved for
equate = np.array(
    [sympy.Eq(ellipse, test_functions[0]), sympy.Eq(ellipse, test_functions[1]), sympy.Eq(ellipse, test_functions[2]),
     sympy.Eq(ellipse, test_functions[3]), sympy.Eq(ellipse, test_functions[4])])

# Calculating the intersection points of the ellipse and the testing functions

# Array that holds the bounds of the integral solved symbolically
bounds_symbolic = np.array([])
for i in range(0, 5):
    bounds_symbolic = np.append(bounds_symbolic, sympy.solve(equate[i], x))

# Array of a-values to plug into the bounds of the integral
a_values = np.array(np.linspace(-10, 10, 201))

# Setting b equal to a constant of 1
b = 1

integrand = np.array([])
for j in range(0, 5):
    integrand = np.append(integrand, (ellipse - test_functions[j]))

# New array with a-values substituted into the bounds
bounds_a = bounds_symbolic
# for j in range(0, 5):
#    bounds_a = np.append[:, ]

Thank you!

Answer 1

Numpy arrays are the best choice when working with pure numerical data, for which they can help speed up many types of calculations. Once you start mixing sympy expressions, things can get very messy. You'll also lose all the speed advantages of numpy arrays.

Apart from that, np.append is a very slow operation as it needs to recreate the complete array every time it is executed. When creating a new numpy array, the recommended way it to first create an empty array (e.g. with np.zeros()) already with its final size.

You should also check out Python's list comprehension as it eases the creation of lists. In "pythonic" code, indices are used as little as possible. List comprehension may look a bit weird when you are used to other programming languages, but you quickly get used to them, and from then on you'll certainly prefer them.

In your example code, numpy is useful for the np.linspace command, which creates an array of numbers (again converting them with np.array isn't necessary). And at the end, you might want to convert the substituted values to a numpy array. Note that this won't work when solve would return a different number of solutions for some of the equations, as numpy arrays need an equal size for all its elements. Also note that an explicit conversion from sympy's numerical type to a dtype understood by numpy might be needed. (Sympy often works with higher precision, not caring for the loss of speed.)

Also note that if you assign b = 1, you create a new variable and lose the variable pointing to the sympy symbol. It's recommended to use another name. Just writing b = 1 will not change the value of the symbol. You need subs to substitute symbols with values.

Summarizing, your code could look like this:

import sympy
import numpy as np

a, b, x, y = sympy.symbols("a b x y")

# Equation of the ellipse solved for y
ellipse = sympy.sqrt((b ** 2) * (1 - ((x ** 2) / (a ** 2))))

# Functions to be tested
test_functions = [a * b * x, ((a * b) ** 2) * x, ((a * b) ** 3) * x, ((a * b) ** 4) * x, ((a * b) ** 5) * x]

# Equating ellipse and test_functions so their intersection can be symbolically solved for
# Array that holds the bounds of the integral solved symbolically
bounds_symbolic = [sympy.solve(sympy.Eq(ellipse, fun), x) for fun in test_functions]

# Array of a-values to plug into the bounds of the integral
a_values = np.linspace(-10, 10, 201)

# Setting b equal to a constant of 1
b_val = 1

# New array with a-values substituted into the bounds
bounds_a = [[[bound.subs({a: a_val, b: b_val}) for bound in bounds]
             for bounds in bounds_symbolic]
            for a_val in a_values]
bounds_a = np.array(bounds_a, dtype='float')  # shape: (201, 5, 2)

The values of the resulting array can for example be used for plotting:

import matplotlib.pyplot as plt

for i, (test_func, color) in enumerate(zip(test_functions, plt.cm.Set1.colors)):
    plt.plot(a_values, bounds_a[:, i, 0], color=color, label=test_func)
    plt.plot(a_values, bounds_a[:, i, 1], color=color, alpha=0.5)
plt.legend()
plt.margins(x=0)
plt.xlabel('a')
plt.ylabel('bounds')
plt.show()

Or filled:

for i, (test_func, color) in enumerate(zip(test_functions, plt.cm.Set1.colors)):
    plt.plot(a_values, bounds_a[:, i, :], color=color)
    plt.fill_between(a_values, bounds_a[:, i, 0], bounds_a[:, i, 1], color=color, alpha=0.1)