Getting an Error When Trying to Solve Systems of ODEs

Question

Here is my code:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math

def VectorField(w, t, p):
    """
    Defines the differential equations to find a in terms of eta.

    Arguments:
        w :  vector of the state variables:
                 w = [eta, a]
        t :  time
        p :  vector of the parameters:
                 p = [H_0, m, r, d]
    """
    eta, a = w
    H_0, m, r, d = p

    # Create f = (eta', a'):
    f = [a, H_0 * (m*a^(-1) + r*a^(-2) + d* a^2)^(1/2)]
    return f

# Parameter values
H_0 = 0.000000000000000002184057032
m = 0.315
r = 0.0000926
d = 0.685

# Initial Conditions
a0 = eta_0 = 0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 250

#x axis values
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
initials = [eta_0, a0]
p = [H_0, m, r, d]

# Call the ODE solver.
wsol = odeint(VectorField, initials, t, args=(p,), atol=abserr, rtol=relerr)

However, when I run the code, I can the error message

TypeError: ufunc 'bitwise_xor' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

The code points to my vector f. I am not sure what I did wrong when declaring this vector f. Can anyone tell me what I did incorrectly? How can I get rid of this error? Thanks

Answer 1

This one trips me up sometimes when I go back and forth between languages.

This line, for example:

f = [a, H_0 * (m*a^(-1) + r*a^(-2) + d* a^2)^(1/2)]

should be

f = [a, H_0 * (m*a**(-1) + r*a**(-2) + d* a**2)**(1/2)]

Because ** is the power operator, not ^.