Avoiding hard coding a lot of coupled ODEs in python

Question

First of, I'm sorry if the title is not entirely fitting, I had a hard time finding an appropriate one (which might have also effect my searching efficiency for already asked questions like this :/ ).

The problem is the following. While it is comparably easy to solve coupled ODE's in python with Scipy, I still have to write down my ODE in the form explicitly. For example for a coupled ODE of the form

d/dt(c_0)=a(c_0)+b(c_1) and d/dt(c_1)=c(c_0)

I would set up sth like:

import numpy as np
from scipy.integrate import ode

a=1
b=2
c=3
val=[]

def dC_dt(t, C):
    return [a*C[0]+b*C[1],
            c*C[0]]

c0, t0 = [1.0,0.0], 0
r = ode(dC_dt).set_integrator('zvode', method='bdf',with_jacobian=False)
r.set_initial_value(c0, t0)
t1 = 0.001
dt = 0.000005
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
val.append(r.y)

However, now I have coupled ODE's of the rough form

d/dt(c_{m,n})=a(c_{m,n})+b(c_{m+1,n-1})+k(c_{m-1,n+1})

with c_{0,0}=1 and I have to include orders with m^2+n^2-mn smaller than a max value. For a small max, what I did, is using a dictionary to use a notation with two indices and map it on a 1D list

dict_in={'0,0':0,'-1,0':2,...}

and then I entered the ODE for each order

def dC_dt(t,C):
    return[a*C[dict_in['0,0']]+b*C[dict_in['1,-1']]...

Now I basically have to do that for some 100 coupled equations, which I ofc do not want to hard code, so I was trying to figure out a way, to realize the ODE's with a loop or sth. However I couldn't yet find a way around the fact of having two indices in my coefficients together with the condition of only including orders with m^2+n^2-mn smaller than a max value. As I am running in some deadlines, I figured it is time to ask smarter people for help.

Thanks for reading my question!

So c[m,n] depends on c[m+1,n-1] which in turn is influenced by c[m+2,n-2], and so on for every c[m+k,n-k], and also c[m-k,n+k]. Is there any stopping condition, any largest k? And if yes, what is the boundary condition? Is this all the structure of your problem, does the set of coordinates resp. gridpoints with non-zero c value stay the same over all times? — Lutz Lehmann, Apr 15 '15 at 06:58

score 1 · Answer 1 · 2015-04-15T16:28:24.263

1

I had a similar problem. If you fill you dictionary you can just redeclare the function more times inside the loop. This is a silly example of how it works:

dict_in={'0,0':0,'-1,0':2}

for elem in dict_in:
    def dC_dt(t,C):
        #return[a*C[dict_in['0,0']]+b*C[dict_in['1,-1']]
        return dict_in[elem]

    t, C = 0, 0
    print(dC_dt(t,C))
    #r = ode(dC_dt).set_integrator('zvode', method='bdf',with_jacobian=False)

If you need to use more functions together you can use anonymous functions and store them in memory. Another example:

functions_list = list()
for i in range(4):
    f = lambda n = i:  n
    functions_list.append(f)

for j in range(4):
    print(functions_list[j]())

You can use a list or a generator too. For example you can write down the value on a txt file and read that with the readline function each time.

As pointed in the comments below, if you use lamda functions you should pay attention to references. See also https://docs.python.org/3/faq/programming.html#why-do-lambdas-defined-in-a-loop-with-different-values-all-return-the-same-result

edited Apr 15 '15 at 16:28

answered Apr 14 '15 at 20:02

Thanks for your answer. I am not sure I understood it, though. Wouldn't the constant redeclaring of the function dC_dt within a loop of the elements lead to some sort of "uncoupling" of the different elements? Since I use the ode function on exactly one dC_dt, I do require it to update each of the elements of the dictionary at each step of the integration, since each element also depends on the time development of other elements (which are not necessarily its neighbors in the list, which adds a further complication). – SevD Apr 14 '15 at 23:08
Actually I am not sure I understood your problem. I thought you would use each function only one time. In that case, you can declare only one function and change its body each time, like I did in the first example. If you need instead to use each function more than one time use the other way and store them in memory like in the second example. – Apr 15 '15 at 13:14
In your last comment you said you need to change the elements of the dictionary. Why do you need to update them if they are boundary conditions? – Apr 15 '15 at 13:38
1

The second example above is not correct, because of http://stackoverflow.com/questions/938429/scope-of-python-lambda-functions-and-their-parameters – pv. Apr 15 '15 at 15:16
+1. You are right, that example works but is misleading. My point was to show only how to create functions with a loop and I didn't see the issue. I'll change that. – Apr 15 '15 at 16:24
Not "references", but **late binding**. – Karl Knechtel Aug 19 '22 at 03:55

Avoiding hard coding a lot of coupled ODEs in python

1 Answers1