Do calculation for array of values - Python
Question: I am trying to calculate and plot the amplitude `amp` of the solution of a
time dependent `t` differential equation of motion (see `rhs6`) as a function
of `wd` for multiple values of the force coefficient `f` found in `f_array`.

So far I have got the code working for plotting `amp` against `wd` for a
single value of `f`. The result is a resonance peak:

![Resonance peak](http://i.stack.imgur.com/EkhUN.png)

The code for plotting `amp` against `wd` for one value of `f` (which produced
the above image) is given below.

    
    
    from pylab import *
    from scipy.integrate import odeint
    from numpy import *
    import math
    
    #Parameters.
    k = 2.0
    m = 1.0
    w0 = (k/m)**(1/2)
    alpha = 0.2
    l = alpha/(2*m)
    f = 1.0
    wd = w0 + 0.025
    beta = 0.2
    t_fixed = 2.0
    
    #Arrays.
    t = linspace(0.0, 400.0, 400.0)
    wd_array = linspace(w0-1.0, w0+1, 400.0)
    f_array = linspace(10.0, 100.0, 3.0)
    
    #Time step.
    init_x = 0.0
    init_v = 0.0
    dx = 15.0
    dv = 0.0
    init_cond = [init_x,init_v]
    init_cond2 = [init_x + dx,init_v + dv]
    
    def rhs6(c,t,wd):
        c0dot = c[1]
        c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t)
        return [c0dot, c1dot]
    
    amp_array=[]
    for wd in wd_array:
        res = odeint(rhs6, init_cond, t, args=(wd,))
        amp = max(res[:,0]) 
        amp_array.append(amp)
    
    plot(wd_array, amp_array)
    xlabel('Driving frequency, wd')
    ylabel('Ampltiude, amp')    
    show()
    

Now I wish to find `amp` against `wd` for multiple values of `f`. I have tried
to do this by basically making a `for` loop statement over the `f_array`.
However, my approach does not work and I get the error:

`setting an element with a sequence`.

As it's good to show an attempt, below is mine.

    
    
    from pylab import *
    from scipy.integrate import odeint
    from numpy import *
    import math
    
    #Parameters.
    k = 2.0
    m = 1.0
    w0 = (k/m)**(1/2)
    alpha = 0.2
    l = alpha/(2*m)
    f = 1.0
    wd = w0 + 0.025
    beta = 0.2
    t_fixed = 2.0
    
    #Arrays.
    t = linspace(0.0, 200.0, 200.0)
    wd_array = linspace(w0-1.0, w0+1, 200.0)
    f_array = linspace(10.0, 200.0, 3.0)
    
    #Time step.
    init_x = 0.0
    init_v = 0.0
    dx = 15.0
    dv = 0.0
    init_cond = [init_x,init_v]
    init_cond2 = [init_x + dx,init_v + dv]
    
    def rhs6(c,t,wd,f):
        c0dot = c[1]
        c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t)
        return [c0dot, c1dot]
    
    full_array = zeros(len(f_array))    
    for index,force in enumerate(f_array):
        amp_list = []
        for wd in wd_array:
            res = odeint(rhs6, init_cond, t, args=(wd,force))
            amp = max(res[:,0]) 
            amp_list.append(amp)
            print(res)
        amp_array = array(amp_list)
        full_array[index] = amp_array
    
    for f in full_array:
        plot(wd, amp)
    
    show()
    

Any ideas?


Answer: Your issue is that `full_array` is a numpy array and you're trying to set an
element inside it as a list, hence the `setting an element with a sequence`.

To solve this you can instead create a two-dimensional numpy array and then
set each of the rows to be `amp_array` as below

    
    
    from pylab import *
    from scipy.integrate import odeint
    from numpy import *
    import math
    
    #Parameters.
    k = 2.0
    m = 1.0
    w0 = (k/m)**(1/2)
    alpha = 0.2
    l = alpha/(2*m)
    f = 1.0
    wd = w0 + 0.025
    beta = 0.2
    t_fixed = 2.0
    
    #Arrays.
    t = linspace(0.0, 200.0, 200.0)
    wd_array = linspace(w0-1.0, w0+1, 200.0)
    f_array = linspace(10.0, 200.0, 3.0)
    
    #Time step.
    init_x = 0.0
    init_v = 0.0
    dx = 15.0
    dv = 0.0
    init_cond = [init_x,init_v]
    init_cond2 = [init_x + dx,init_v + dv]
    
    def rhs6(c,t,wd,f):
        c0dot = c[1]
        c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t)
        return [c0dot, c1dot]
    
    full_array = zeros((len(f_array),len(wd_array)))
    for index,force in enumerate(f_array):
        amp_list = []
        for wd in wd_array:
            res = odeint(rhs6, init_cond, t, args=(wd,force))
            amp = max(res[:,0]) 
            amp_list.append(amp)
            # print(res)
        amp_array = array(amp_list)
        full_array[index,:] = amp_array
    
    for f in full_array:
        plot(wd_array, f)
    
    show()
    

Alternatively you could follow
[gboffi's](http://stackoverflow.com/a/26013620/3005188) advice and use either
[`np.vstack`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.vstack.html)
or
[`np.hstack`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.hstack.html)
to create your `full_array`.

Your `for` loop at the end to plot the arrays was also incorrect and so I've
edited that. The graph looks like ![Plot](http://i.stack.imgur.com/6tQdf.png)

