A l oop для записи уравнений, которые будут использоваться с odeint - PullRequest
Новая
       75

A l oop для записи уравнений, которые будут использоваться с odeint

0 голосов
/ 14 апреля 2020

У меня есть начальная проблема, которую нужно решить; дифференциальные уравнения получены из словаря, который выглядит следующим образом: 

eqs = {'a': array([-1.,  2.,  4.,  0., ...]),
 'b': array([  1., -10.,   0.,   0., ...]),
 'c': array([ 0.,  3., -4.,  0., ...]),
 'd': array([ 0.,  5.,  0., -0., ...]),
...}

Дифференциальное уравнение da/dt задано как -1*[a]+2*[b]+4*[c]+0*[d] ....

Используя приведенный выше словарь, я пишу функция dXdt as: 

def dXdt (X, t):
    sys_a, sys_b, sys_c, sys_d,... = eqs['a'], eqs['b'], eqs['c'], eqs['d'],...
    dadt = sys_a[0]*X[0]+sys_a[1]*X[1]+sys_a[2]*X[2]+sys_a[3]*X[3]+...
    dbdt = sys_b[0]*X[0]+sys_b[1]*X[1]+sys_b[2]*X[2]+sys_b[3]*X[3]+...
    dcdt = sys_c[0]*X[0]+sys_c[1]*X[1]+sys_c[2]*X[2]+sys_c[3]*X[3]+...
    dddt = sys_d[0]*X[0]+sys_d[1]*X[1]+sys_d[2]*X[2]+sys_d[3]*X[3]+...
    ...
    return [dadt, dbdt, dcdt, dddt, ...]

Начальные условия: 

X0 = [1, 0, 0, 0, ...]

, а решение имеет вид: 

X = integrate.odeint(dXdt, X0, np.linspace(0,10,11))

Это хорошо работает для небольшая система, где я могу написать уравнения вручную. Однако у меня есть система, которая имеет ~ 150 дифференциальных уравнений, и мне нужно автоматизировать способ, которым я пишу dXdt, чтобы использовать with scipy.integrate.odeint, учитывая словарь eqs. Есть ли способ сделать это?

1 Ответ

1 голос
/ 14 апреля 2020

Каждый раз, когда что-то следует простому линейному шаблону, вы можете использовать итерацию или понимание для express этого. Если у вас есть несколько таких шаблонов, вы можете просто вложить их. Итак, это: 

    sys_a, sys_b, sys_c, sys_d,... = eqs['a'], eqs['b'], eqs['c'], eqs['d'],...
    dadt = sys_a[0]*X[0]+sys_a[1]*X[1]+sys_a[2]*X[2]+sys_a[3]*X[3]+...
    dbdt = sys_b[0]*X[0]+sys_b[1]*X[1]+sys_b[2]*X[2]+sys_b[3]*X[3]+...
    dcdt = sys_c[0]*X[0]+sys_c[1]*X[1]+sys_c[2]*X[2]+sys_c[3]*X[3]+...
    dddt = sys_d[0]*X[0]+sys_d[1]*X[1]+sys_d[2]*X[2]+sys_d[3]*X[3]+...
    ...
    [dadt, dbdt, dcdt, dddt, ...]

можно выразить просто как: 

[sum(eqs[char][i] * X[i] for i in range(len(X))) for char in eqs.keys()]
...