current community

your communities

more stack exchange communities

company blog
Stack Exchange Inbox Reputation and Badges
tour help
Join the Stack Overflow Community
Stack Overflow is a community of 6.4 million programmers, just like you, helping each other.
Join them; it only takes a minute:
Sign up
up vote 0 down vote favorite

I'm quite new to Python and Numpy, so I apologize if I'm missing something obvious here.

I have a function that solves a system of 2 differential equations :

import numpy as np
import numpy.linalg as la

def solve_ode(x0, a0, beta, t):
    At = np.array([[0.23*t, (-10**5)*t], [0, -beta*t]], dtype=np.float32)

    # get eigenvalues and eigenvectors
    evals, V = la.eig(At)
    Vi = la.inv(V)

    # get e^At coeff
    eAt = V @ np.exp(evals) @ Vi

    xt = eAt*x0
    return xt

However, running it with this code :

import matplotlib.pyplot as plt

# initial values
x0 = 10**6
a0 = 2.5
beta = 0.05

t = np.linspace(0, 3600, 360)

plt.semilogy(t, solve_ode(x0, a0, beta, t))

... throws this error :

ValueError: setting an array element with a sequence.

At this line :

At = np.array([[0.23*t, (-10**5)*t], [0, -beta*t]], dtype=np.float32)

Note that t and beta are supposed to be floats. I think Python might not be able to infer this but I don't know how I could do this...

Thx in advance for your help.

share|improve this question
    
If t is a scalar At becomes a 2x2 array. Now the problem is that Numpy doesn't know what to do when you make t a list of numbers. It can't fill the cells with arrays. That why you get the Value Error. Since I'm not a Physics/Maths major I'm not sure how the appropriate solution would look like. Is the idea to run solve_ode() for each value of t? Or should At grow in a specific manner? – Peter Smit yesterday
up vote 1 down vote accepted

You are supplying t as a numpy array of shape 360 from linspace and not simply a float. The resulting At numpy array you are trying to create is then ill formed as all columns must be the same length. In python there is an important difference between lists and numpy arrays. For example, you could do what you have here as a list of lists, e.g.

At = [[0.23*t, (-10**5)*t], [0, -beta*t]]

with dimensions [[360 x 360] x [1 x 360]].

Alternatively, if all elements of At are the length of t the array would work,

At = np.array([[0.23*t, (-10**5)*t], [t, -beta*t]], dtype=np.float32)

with shape [2, 2, 360].

share|improve this answer
    
Thank you very much! For some stupid reason, I assumed scipy would iterate over t and pass to my function every single value contained in the linspace. I understand now, thanks a lot! – Scrashdown yesterday
up vote 1 down vote

When you give a list or a list of lists, or in this case, a list of list of listss, all of them should have the same length, so that numpy can automatically infer the dimensions (shape) of the resulting matrix.

In your example, it's all correctly put, except the part you put 0 as a column I guess. Not sure what to call it though, cause your expected output is a cube I suppose.

You can fix it by giving the correct number of zeros as bellow:

At = np.array([[0.23*t, (-10**5)*t], [np.zeros(len(t)), -beta*t]], dtype=np.float32)

But check the .shape of the resulting array, and make sure it's what you want.

share|improve this answer
up vote 0 down vote

As others note the problem is the 0 in the inner list. It doesn't match the 360 length arrays generated by the other expressions. np.array can make an object dtype array from that (2x2), but can't make a float one.

At = np.array([[0.23*t, (-10**5)*t], [0*t, -beta*t]])

produces a (2,2,360) array. But I suspect the rest of that function is built around the assumption that At is (2,2) - a 2d square array with eig, inv etc.

What is the return xt supposed to be?

Does this work?

S = np.array([solve_ode(x0, a0, beta, i) for i in t])

giving a 1d array with the same number of values as in t?

I'm not suggesting this is the fastest way of solving the problem, but it's the simplest, especially if you are only generating 360 values.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.