issue with converting rate equation to python code

Question

I am trying to convert these rate equations to python code, I have made I lot of research but can't seem to get any clear path to follow to achieve this, please any help will be appreciated

This is a newly updated code....i wrote using the quide from Tom10.....please what do you think?

    import numpy as np
# import numpy as sum  # not necessary, just for convenience, and replaces the builtin

# set N_core value
N_CORE = 0

# set the initial conditions appropriately (you need to set these correctly)
N = np.ones(8)
r = np.ones((8, 8))
dN = np.zeros(8) # the value here is not important for your equations

# set constant for equation 1
R_P1abs37 = 20
F_P1 = 20
R_P1abs47 = 40
W_3317 = 1.0

# set constant for equation 2
W_6142 = 90
W_5362 = 80

# Set you constants appropriately for equation 3
R_P2abs35 = 30
F_P2 = 40
R_L2se34 = 50
F_L2 = 90

# equation 4 constants
W_2214 = 20

#equation 5 constants
R_P1abs13 = 30
R_L2se32 = 20
F_L1 = 10

# equation 1 formular
dN[7] =sum(r[7,:]*N[7]) + (R_P1abs37*F_P1) + (R_P1abs47*F_P1) + (W_3317*N[3]**2)

# equation 2 formular
dN[6] = (r[7,6]*N[7]) - sum(r[6,:]*N[6]) - (W_6142*N[6]*N[1]) + (W_5362*N[5]*N[3])

#equation 3 formular
dN[5] = sum(r[:,5]*N) - sum(r[5,:]*N[5]) + R_P2abs35*F_P2 - R_L2se34*F_L2 - W_5362*N[5]*N[3]

# equation 4 formular
dN[4] = sum(r[:,4]*N) - sum(r[4,:]*N[4]) - (R_P1abs47*F_P1) + (R_L2se34*F_L2) + (W_2214*N[2]**2)+ (W_6142*N[6]*N[1])

#equation 5 formular
dN[3] = sum(r[:,3]*N) - sum(r[3,:]*N[3]) + (R_P1abs13*F_P1) - (R_P1abs37*F_P1) - (R_P2abs35*F_P2)
-(R_L2se32*F_L1) - ((2*W_3317)*N[3]**2) - (W_5362*N[5]*N[3])

#equation 6 formular
dN[2] = sum(r[:,2]*N) - (r[2,1]*N[2]) + (R_L2se32*F_L1) - ((2*W_2214)*N[2]**2) + (W_6142*N[6]*N[1])+(W_5362*N[5]*N[3])


#equation 7 formular
dN[1] = sum(r[:,1] * N) - (R_P1abs13*F_P1) + (W_2214*N[2]**2) + (W_3317+N[3]**2) - (W_6142+N[6]*N[1])

#equation for N CORE
N_CORE = sum(dN)


print(N_CORE)

to be able to pass in some constants and value and compute them with the equations in python — kayode olayiwola, Oct 04 '20 at 17:25
To explain this from the ground up is very difficult. Do you know what a variable is in math and programming? Do you know that an equation is in math and an assignment statement is in programming and do you know the difference? Do you know what the sigma symbol means, and all the subscripts on the variables? If you give a start of what you've done so far, we will have some idea of what you know and we will know where to start an answer. — tom10, Oct 04 '20 at 18:32
@tom10 i have a good knowledge of all what you asked...and i am an intermediate python programmer...the issue i am having current is interpreting those equations.......like......does that summation symbol mean i have to execution the formulae in front of it in a loop of n-times?.....are components like w(subscipt)3317 and Rp1 treated as a variable? — kayode olayiwola, Oct 04 '20 at 18:46

tom10 · Accepted Answer · 2020-10-08T15:36:20.907

0

Here is list of relevant issues based on your question and comments:

Usually if the summation is over i, then everything without an i subscript is constant for that sum. (Mathematically these constant terms can just be brought out of the sum; so the first equation is a bit odd where the N_7 could be moved out of the sum but I think they're keeping it in to show the symmetry with the other equations which all have an r*N term).
The capitol sigma symbol (Σ) means you need to do a sum, which you can do in a loop, but both Python list and numpy have a sum function. Numpy has the additional advantage that multiplication is interpreted as multiplication of the individual elements, making the expression easier. So for a[0]*[b0] + a[1]*b[1] + a[2]*b[2] and numpy arrays is simply sum(a*b) and for Python lists it's sum([a[i]*b[i] for in range(len(a))]

Therefore using numpy, the setup and your third equation would look like:

import numpy as np
import numpy.sum as sum  # not necessary, just for convenience, and replaces the builtin

# set the initial conditions appropriately (you need to set these correctly)
N = np.ones(7, dtype=np.float)
# r seems to be a coupling matrix, and should be set according to your system
r = np.ones((7, 7), dtype = np.float) 
# the values for dN are not important for your equations because dN only appears on the left side of the equations, so we just make a place to store the results
dN = np.zeros(7, dtype=np.float) 

# Set you constants appropriate.y
R_P2abs35 = 1.0
F_P2 = 1.0
R_L2se34 = 1.0
F_L2 = 1.0
W_5362 = 1.0

dN[5] = sum(r[:,5]*N) - sum(r[5,:]*N[5]) + R_P2abs35*F_P2 - R_L2se34*F_L2 - W_5362*N[5]*N[3]

Note that although the expressions in the sums look similar, the first is essentially a dot product between two vectors and the second is a scalar times a vector so N[5] could be taken out of the sum (but I left it there to match the equation).

Final note: I see you're new to S.O. so I thought it would be helpful if I answered this question for you. In the future, please show some attempt at the code -- it really helps a lot.

edited Oct 08 '20 at 15:36

answered Oct 04 '20 at 19:30

tom10

67,082
10
127
137

wooow.....thank you so much......these clarified a lot of things for me......I really appreciate this......I will start the development now.....I will let you know how it turns out.....thank you tom – kayode olayiwola Oct 04 '20 at 20:33
Great. In the end, I hope the question doesn't get closed because I think it's a good topic for a question. Just for general curiosity, I'm interested to know what the equation represents, or at least the field of study, if it's reasonable to state. – tom10 Oct 04 '20 at 20:44
thank you.....can you check the updated question?....i have written the code using your guide......what do you think? – kayode olayiwola Oct 08 '20 at 12:52
@kayodeolayiwola: I spot-checked a few of your equations and they look good to me. The real way to test these is to have an idea of what they should produce in a few situations and make sure they work, and tests similar to that. One mistake I made in my answer is that I should have created the initial arrays as floats not ints, and I've fixed that now in my answer. It's probably also good to set the constants to floats too (eg, `30.0` instead of `30`). Also, in your equations, you'll need to set the values of *r* and the initial values of *N* to match your system. – tom10 Oct 08 '20 at 13:50
ooh...thank you....i have changed it....but what do you mean by setting initial value of r and N?.....do you mean i need to populate the numpy arrays(r and N) with some values before the equations?.....thank you – kayode olayiwola Oct 08 '20 at 14:33
@kayodeolayiwola: Yes, both *r* and *N* need to be set, although they are slightly different. Setting *N* will set the **initial conditions** of the simulation. Most dynamical systems behave in a way that depends on two things: 1) the details of the dynamics of the system, and this is specified by the equations; and, 2) the initial conditions, that is, how it starts off. Think of a pendulum, how far you pull it up before releasing, that's the initial values of *N*, but this is much more complex. – tom10 Oct 08 '20 at 15:00
@kayodeolayiwola: *r* needs to be set too. From looking at the equation, *r* is basically the coupling between the different *N* (although coupling is also added explicilty in a few terms, link N[3]*N[5]). That is, without *r* (and the explicit coupling) these would just be 7 independent equations, but *r* makes each *N* depend on all the others. Eg, once *r* is set, you could see this by starting off with N=(0,0,0,1,0,0,0,0), and eventually all the N will be doing something, but without *r*, only N_4 would ever do anything interesting at all (expect for growing with the constant terms). – tom10 Oct 08 '20 at 15:06
i think i understand this......just to be sure...please can you give me an example using random values on setting the initial condition of N and R?.....thank @tom10 – kayode olayiwola Oct 08 '20 at 15:18
@kayodeolayiwola: Doesn't where you got the equation say something about these, especially *r*? In my answer, the `ones` do this. Anyway, *N* is like `N = np.array([0.0, 3.4, 1.2, 4.5, 0.1, 0.8, 100.1, 4.7]`) (where the first value should be zero so it won't effect the sum). *r* is much harder since it's very big (7x7=49 terms) and reflects the coupling of the system, and would be very tedious to write out. If you're lucky, it's something like the identity with a few side terms, so `r = np.identity(8)` and then a few things like, `r[1,3] = 0.2`, `r[5,2] = 0.007`, etc. – tom10 Oct 08 '20 at 15:46
@kayodeolayiwola: Also, I noticed an error in your equations, the final line should be `N_CORE = sum(N)`, not `dN`. (And, btw, I have to get back to my work now.) – tom10 Oct 08 '20 at 15:48
@kayodeolayiwola: One small error in previous phrasing: *r* seems to be a linear rate (probably why it's called *r*), so it shouldn't be the identity, but maybe a constant times the identity corresponding to the rate, or different terms but mostly along the diagonal of the matrix (whereas the other terms with Ns, like, N[5]*N[3] are non-linear coupling terms). The reason I say prob mostly along the diagonal is that for many systems, how "a thing" (ie, an N) moves depends on how **it** (diagonal) is already moving, and maybe a bit on how other things (off-diagonal) move too. If you're lucky. – tom10 Oct 08 '20 at 15:59

issue with converting rate equation to python code

1 Answers1