Solve ODE for Parameter Arrays (Python)

* I know this question is pretty simple, but I would like to know the best way to set up such a for loop in Python.

I already wrote a program to calculate and construct a solution to a 2nd order differential equation (this code is shown below).

I would like to know the best techniques for iterating this calculation for an array of f parameters (hence f_array

). That is, so the graph shows the 20

data sets related to solutions as a function of t

, each with a different value f

.

Cheers for any ideas.

from pylab import *
from scipy.integrate import odeint

#Arrays.
tmax = 100
t = linspace(0, tmax, 4000)
fmax = 100
f_array = linspace(0.0, fmax, 20)

#Parameters
l = 2.5
w0 = 0.75
f = 5.0
gamma = w0 + 0.05
m = 1.0
alpha = 0.15
beta = 2.5

def rhs(c,t):
    c0dot = c[1]
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((gamma)*t)-alpha*c[0] - beta*c[0]*c[0]*c[0]
    return [c0dot, c1dot]

init_x = 15.0
init_v = 0.0
init_cond = [init_x,init_v]
ces = odeint(rhs, init_cond, t)

s_no = 1
subplot(s_no,1,1)
xlabel("Time, t")
ylabel("Position, x")
grid('on')
plot(t,ces[:,0],'-b')
title("Position x vs. time t for a Duffing oscillator.")
show()

      

Here is a graph showing the solution to this equation with respect to a single value f

for an array of values t

. I would like to quickly repeat this plot for an array of values f

.

http://i61.tinypic.com/28bgyzs.png

+3


source to share


1 answer


Here's one approach:

Modify rhs

to accept a third argument, parameter f

. The definition rhs

must begin

def rhs(c, t, f):
    ...

      

Iterating over f_array

with a loop for

. In the loop, call odeint

with an argument args

to odeint

give the value f

as the third argument rhs

. Save the results of each call odeint

in a list. Basically replace

ces = odeint(rhs, init_cond, t)    

      



from

solutions = []
for f in f_array:
    ces = odeint(rhs, init_cond, t, args=(f,))
    solutions.append(ces)

      

For each value f

in f_array

, you now have a solution listed solutions

.

To build them, you can place your call plot

in another loop for

:

for ces in solutions:
    plot(t, ces[:, 0], 'b-')

      

+4


source







All Articles