CODEWITHSUNDEEP
pythonnumpymathscipynumerical-integration
MD'
MD'

Reputation: 105

implement an integration math equation using odeint in Python

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: Original equation

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:

Validation image

and here is the numpy output:

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?

Upvotes: 3

Views: 386

Answers (1)

Pierre de Buyl
Pierre de Buyl

Reputation: 7293

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, '-')

Upvotes: 1

Related Questions