-3

My math question is here - https://math.stackexchange.com/questions/2063507/solving-this-integral-involving-ei-function

This relates Population dynamics with t = time, N_t = population at time t, r = rate of growth and K = cqrrying capacity.

My code is attached below.

Python is unable to calculate the value of N_next because it's inside the exponential integral function scipy.special.expi(). How can I circumvent this? It's now saying "Scipy has no attribute special" but according to this - https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.expi.html - It should be.

import math
import scipy
import numpy as np
import matplotlib
import matplotlib.pyplot as plt

t_f =100
N_0 = 10
t = []
N_t = [N_0,]
r = 2.5
K = 1000

for i in range(0,100):
    scipy.special.expi(r*N_next/K) = i*math.exp(r) + scipy.special.expi(r/K * N_t[i])
    N_t.append(N_next)
    t.append(i)

plt.plot(t,N_t)
plt.show
Community
  • 1
  • 1
Adeetya
  • 101
  • 9
  • 4
    Your code makes no sense. Are you leaving anything out? Is there a function called `ei` you've defined somewhere? But you are trying to assign to the result of a call, but that isn't possible. Which is why you are getting an error. – juanpa.arrivillaga Dec 21 '16 at 19:20
  • Sorry. I missed a part. 1 sec – Adeetya Dec 21 '16 at 19:20
  • 3
    Also, `N_next` isn't defined. – juanpa.arrivillaga Dec 21 '16 at 19:21
  • 1
    also note that `^` is `xor`; use `math.exp` if you mean that... – hiro protagonist Dec 21 '16 at 19:23
  • N_next is a variable defined inside the loop. I made other programs with this format, and they worked perfectly. – Adeetya Dec 21 '16 at 19:23
  • 3
    `N_next` is not defined anywhere in your code. If your other programs had this format, in Python, I guarantee you they do not work. – juanpa.arrivillaga Dec 21 '16 at 19:24
  • Calculate N_next. – Adeetya Dec 21 '16 at 19:24
  • Thanks @hiroprotagonist - I made a mistake there too I guess – Adeetya Dec 21 '16 at 19:25
  • 1
    `scipy.special.expi(r*N_next/K) = i*math.exp(r) + scipy.special.expi(r/K * N_t[i])` still does not make any sense whatsoever - You seem to be trying to throw a mathematical function at Python and hope it will magically solve it for you – UnholySheep Dec 21 '16 at 19:34
  • 3
    `=` in Python is not for declaring equations but for **assigning** the result of evaluating the **expression on the right-hand side** to the **reference on the left-hand side**. A function call is not a valid reference. (Variable names and indexing into existing variables would be references that you can assign to.) – das-g Dec 21 '16 at 19:37
  • The RHS is an expression where python takes values from the N_t list and the value of i to get a constant, and LhS is where the variable is, N_next – Adeetya Dec 21 '16 at 19:37
  • @das-g So how do I get python to solve this equation? – Adeetya Dec 21 '16 at 19:39
  • 1
    Either solve it yourself (analytically) so that you can give Python just an expression to evaluate, or use a library for solving equations. – das-g Dec 21 '16 at 19:41
  • For the latter see [How to solve a pair of nonlinear equations using Python?](http://stackoverflow.com/questions/8739227/how-to-solve-a-pair-of-nonlinear-equations-using-python) – das-g Dec 21 '16 at 19:44
  • But I don't have two equations or two variables, it's just one variable - N_next – Adeetya Dec 21 '16 at 19:47
  • It might help if you were to show us what mathematical problem you want to solve or explore. Put it in the question. – Bill Bell Dec 21 '16 at 19:52
  • @BillBell - added the link - please have a look – Adeetya Dec 21 '16 at 19:57
  • @BillBell - Will this be easier in a different programming language? – Adeetya Dec 21 '16 at 20:18
  • One _huge_ advantage programming languages have over pen-and-ink mathematical notation is that you don't need to clutter your algorithm with clarifications like "`N_t` = population at time `t`" --- you can just give the variable a _meaningful_ name and initial value, like `pop = [10]`. Also, there's no need to "subscript" variables TeX-style --- array indexing _is_ your subscript: `pop[0]` is identical to `N_0`, with the added advantage of not being wrong if the first value in `pop` ever changes. – Kevin J. Chase Dec 21 '16 at 20:20
  • Frankly, once you become accustomed to Python, you won't find anything much easier. – Bill Bell Dec 21 '16 at 20:30

3 Answers3

0

That error is the Python interpretr fancy way of telling: You cannot assign an expression result to a function call

for i in range(0,100):
    # function call assignment expression
    ei(r*N_next/K) = i*e^r + ei(r/K * N_t[i])

Maybe you mixed the line in your example? In that case please edit your post so we can help you.

yorodm
  • 4,359
  • 24
  • 32
  • This isn't really an answer, it's more of a comment – UnholySheep Dec 21 '16 at 19:35
  • But I need to do that. It's something like ei(m*x) = C.. where m and c are both constants after calculations. And I need to find the value of x. Is it not possible to do that? I have made a few edits regarding scipy.special.expi = ei(in my code). Please check – Adeetya Dec 21 '16 at 19:36
  • 3
    As I said in the answer `ei(r*N_next/K) = i*e^r + ei(r/K * N_t[i])` that line is the problem, you **cannot** assign an expression result to a function call, meaning that the thing on the left side of `=` cannot be a function call. – yorodm Dec 21 '16 at 19:44
0

The answer in https://math.stackexchange.com/questions/2063507/solving-this-integral-involving-ei-function provides a formula that you want to evaluate for a sequence of values of t. In building a script such as this I suggest that you take a bottom-up approach. Here's what I have written so far.

from numpy import exp
import matplotlib.pyplot as plt

K = 1000
r = 2.5

N_0 = 10
N_t = 50

def Ei(x):
    return x

def Ei_inv(x):
    return 1/x

def N(t):
    return K*Ei_inv(exp(r)*t+Ei(r*N_0/K))/r

t = [_ for _ in range(40)]
N = [N(_) for _ in t]

plt.plot(t,N)
plt.show()
  • I knew what some of the constants are; I therefore put them in the code.
  • I didn't know what Ei or its inverse would be, or where I might find these functions; as an interim measure I wrote 'dummy' functions to take their places.
  • Now I could write the function that I actually wanted to evaluate as another Python function with a single parameter, in terms of the dummy functions. At that point I made a trial execution of what I had using print (N(50)) to verify that the various bits of code were behaving politely together.
  • I added two lines, one to define the x, or time, values for a plot, and the second the y, or N values.

Your turn, replace the dummy functions with whatever is needed to calculate Ei and its inverse.

Community
  • 1
  • 1
Bill Bell
  • 21,021
  • 5
  • 43
  • 58
0

So, I figured it out. Well,kinda. The code is now very computationally intensive. But it works atleast. Here it is below.

import math
from scipy import *
from scipy.special import expi
import numpy as np

t_f =100
N_0 = 10
t = [0,]
n = [N_0,]
r = 2.5
K = 100
N_t = [N_0,]

for i in range(0,100):
    a= i*math.exp(r) + expi(r/K * N_0)
    n.append(a)
    t.append(i+1)

g = np.arange(0,1000,0.0001)
d = []
e = []

for i in g:
    x = expi(i)
    d.append(x)
    e.append(i)    

for j in range(1,99):
    b = n[j]
    c = []

    for i in g:
        c.append(b)

    t = np.subtract(d,c)
    for i in range(len(t) - 1):
        if t[i] == 0. or t[i] * t[i + 1] < 0.:
            y = e[i]
            print(y)

    N_next = math.floor(y*K/r)
    N_t.append(N_next)

print(N_t)
plot(t,N_t)
Adeetya
  • 101
  • 9