implement an integration math equation using odeint in Python

Question

I am trying to solve the following equation in python using the scipy.odeint function.

equation 1

Currently I am able to implement this form of the equation

equation 2

in python using the following script:

def dY(y1, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    return (a/dC)*(yin-y1)+b1*dC

x = np.linspace(0,20,1000)
y0 = 0
res = odeint(dY, y0, x)
plt.plot(t,res, '-')
plt.show()

My problem with the first equation is 'i'. I don't know how to integrate the equation and still be able to provide the current and previous 'y'(yi-1 and yi) values. 'i' is simply a sequence number that is within a range of 0..100.

Edit 1:

The original equation is:

Which I rewrote using y,x,a,b and C

Edit2: I edited Pierre de Buyl' code and changed the N value. Luckily I have a validation table to validate the outcome against. Unfortunately, the results are not equal.

Here is my validation table:

and here is the numpy output:

Used code:

def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 3
    dC = C/N
    b1 = 0.01
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    return (a/dC)*(y_diff)+b1*dC

x = np.linspace(0,20,11)
y0 = np.zeros(3)
res = odeint(dY, y0, x)
plt.plot(x,res, '-')

as you can see the values are different by an offset of 0.02..

Am I missing something that results in this offset?

That notation could mean that you're using Euler forward integration for a single ODE for y OR integrating multiple coupled ODEs. Which one is it? — duffymo
– duffymo, Commented Sep 14, 2017 at 17:58
@duffymo they are using the Euler forward integration for solving y — MD'
– MD', Commented Sep 14, 2017 at 18:51
@MD' is this differential equation correct? Could you please give a reference from where you got it? — Saullo G. P. Castro
– Saullo G. P. Castro, Commented Sep 15, 2017 at 0:36
@SaulloCastro I added the original equation to the original post. The equation is calculating the concentration of a reactant in a pipe reactor. — MD'
– MD', Commented Sep 15, 2017 at 6:51
Is this your equation? cee.mtu.edu/~reh/courses/ce251/251_notes_dir/node3.html — duffymo
– duffymo, Commented Sep 15, 2017 at 12:06

Pierre de Buyl · Accepted Answer · 2017-09-15 08:01:28Z

1

The equation is a "coupled" ordinary differential equation (see "System of ODEs" on Wikipedia.

The variable is a vector containing y[0], y[1], etc. To solve the ODE you must feed a vector as the initial condition and the function dY must return a vector as well.

I have modified your code to achieve this result:

def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    return (a/dC)*y_diff+b1*dC

I have written the part y[i-1] - y[i] as a NumPy vector operation and special cased the coordinate y[0] (that is the y1 in your notation but arrays start at 0 in Python).

The solution, using an initial value of 0 for all yi is

x = np.linspace(0,20,1000)
y0 = np.zeros(4)
res = odeint(dY, y0, x)
plt.plot(x,res, '-')

answered Sep 15, 2017 at 8:01

Pierre de Buyl

7,3232 gold badges18 silver badges23 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

MD' Over a year ago

thank you for the implementation. However, I executed the code with a little modification to y0 and N such that the output matches the results of my validation table. See Edit2 where I elaborate on what I have.

Collectives™ on Stack Overflow

implement an integration math equation using odeint in Python

1 Answer 1

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related