python求积分的几种原理,在Python中使用odeint实现积分数学方程

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

Currently I am able to implement this form of the equation

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?

解决方案

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