Writing functions from dictionaries to be used with solve_ivp

Question

I am trying to solve a large system of differential equations using solve_ivp.

from scipy import integrate

def dXdt(X,t):
    return np.array([dadt(X,t), dbdt(X,t), dcdt(X,t), dddt(X,t])

sol = integrate.solve_ivp(dXdt, (0,100), initial_value_array, t_eval)

The dadt(X,t), dbdt(X,t), dcdt(X,t), dddt(X,t] is a system of differential equations that I need to obtain from the following dictionaries:

da_dict = {'a': -1.0, 'b': 2.0, 'c': 4.0}
db_dict = {'b': -10.0, 'a': 1.0}
dc_dict = {'c': -4.0, 'b': 3.0}
dd_dict = {'b': 5.0}

as follows:

def dadt(X,t):
    return -1.0*X[0] + 2*X[1] + 3*X[2]

where, X[0], X[1], X[2], x[3] are represented by 'a', 'b', 'c', 'd' in the dictionaries. Similarly,

def dbdt(X,t):
    return -10*X[1] + 1*X[0]

def dcdt(X,t):
    return -4*X[2] + 3*X[1]

def dddt(X,t):
    return 5*X[1]

I have more than 100 differential equations I need to solve using solve_ivp. How do I write dadt(X,t), dbdt(X,t)....from the dictionaries without actually writing them?

What is the issue, exactly? Have you tried anything, done any research? See [ask], [help/on-topic]. — AMC, Apr 07 '20 at 19:09
Is the large system always linear? Then you might get better results using the matrix exponential as implementation of the exact solution. I'm not sure if there are sparse variants. — Lutz Lehmann, Apr 07 '20 at 19:39
@LutzLehmann, It is not. The system is a combination of first and second order differential equations. — DPdl, Apr 07 '20 at 19:42
Are the expressions encoded in the python source as dictionaries? Then you could with the same effort encode them as sympy expressions and apply the code-generation facilities like autowrap or ufuncify to generate the ODE function. Or use jitcode which does it more automatically. — Lutz Lehmann, Apr 08 '20 at 07:55

score 0 · Answer 1 · answered Apr 07 '20 at 18:41

I'm not sure what to do with the t parameter since you didn't use it in any of the examples. I put placeholder values for X and t in the code below

I wrote a snippet below that could be useful for your example. I put your dictionaries in a list if you wrote all of them already you can add the others yourself

da_dict = {'a': -1.0, 'b': 2.0, 'c': 4.0}
db_dict = {'b': -10.0, 'a': 1.0}
dc_dict = {'c': -4.0, 'b': 3.0}
dd_dict = {'b': 5.0}

# associate the order of the coefficient letter to a position in the array.
order = {
    "a": 0,
    "b": 1,
    "c": 2,
    "d": 3,
    "e": 4
}

def equation_maker(X,t,equation):

    output_eq = 0
    # coeff: a,b,.. coeff_value: -1.0, 2.0,...
    for coeff,coeff_value in equation.items():
        output_eq += X[order.get(coeff)]*coeff_value
    return output_eq


equation_list = [da_dict,db_dict,dc_dict,dd_dict]
# output list of equations that you can convert to numpy arrays or whatever
diff_eq_list = []
# placeholder values for X and t
X = [0,1,2,3,4]
t = 0
for equation in equation_list:
    # assuming you have X and t from somewhere else
    diff_eq_list.append(equation_maker(X,t,equation))

print(diff_eq_list)

I have never used scipy (yet!) but I hope this can help you not having to write all the functions, let me know if this looks useful to you.

Writing functions from dictionaries to be used with solve_ivp

1 Answers1