Odeint python как работает

Численное решение обыкновенных дифференциальных уравнений (ОДУ) в Python

Модуль scipy.integrate имеет две функции ode() и odeint(), которые предназначены для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (т.е. задача Коши).

Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат

Решение одного ОДУ

Допустим надо решить диф. уравнение 1-го порядка

Получилось что-то такое:

Решение системы ОДУ

Пусть теперь мы хотим решить (автономную) систему диф. уравнений 1-го порядка

Integration ( scipy.integrate )#

The scipy.integrate sub-package provides several integration techniques including an ordinary differential equation integrator. An overview of the module is provided by the help command:

General integration ( quad )#

The function quad is provided to integrate a function of one variable between two points. The points can be \(\pm\infty\) ( \(\pm\) inf ) to indicate infinite limits. For example, suppose you wish to integrate a bessel function jv(2.5, x) along the interval \([0, 4.5].\)

This could be computed using quad :

The first argument to quad is a “callable” Python object (i.e., a function, method, or class instance). Notice the use of a lambda- function in this case as the argument. The next two arguments are the limits of integration. The return value is a tuple, with the first element holding the estimated value of the integral and the second element holding an upper bound on the error. Notice, that in this case, the true value of this integral is

is the Fresnel sine integral. Note that the numerically-computed integral is within \(1.04\times10^<-11>\) of the exact result — well below the reported error bound.

If the function to integrate takes additional parameters, they can be provided in the args argument. Suppose that the following integral shall be calculated:

This integral can be evaluated by using the following code:

Infinite inputs are also allowed in quad by using \(\pm\) inf as one of the arguments. For example, suppose that a numerical value for the exponential integral:

is desired (and the fact that this integral can be computed as special.expn(n,x) is forgotten). The functionality of the function special.expn can be replicated by defining a new function vec_expint based on the routine quad :

The function which is integrated can even use the quad argument (though the error bound may underestimate the error due to possible numerical error in the integrand from the use of quad ). The integral in this case is

This last example shows that multiple integration can be handled using repeated calls to quad .

General multiple integration ( dblquad , tplquad , nquad )#

The mechanics for double and triple integration have been wrapped up into the functions dblquad and tplquad . These functions take the function to integrate and four, or six arguments, respectively. The limits of all inner integrals need to be defined as functions.

An example of using double integration to compute several values of \(I_\) is shown below:

As example for non-constant limits consider the integral

This integral can be evaluated using the expression below (Note the use of the non-constant lambda functions for the upper limit of the inner integral):

For n-fold integration, scipy provides the function nquad . The integration bounds are an iterable object: either a list of constant bounds, or a list of functions for the non-constant integration bounds. The order of integration (and therefore the bounds) is from the innermost integral to the outermost one.

The integral from above

can be calculated as

Note that the order of arguments for f must match the order of the integration bounds; i.e., the inner integral with respect to \(t\) is on the interval \([1, \infty]\) and the outer integral with respect to \(x\) is on the interval \([0, \infty]\) .

Non-constant integration bounds can be treated in a similar manner; the example from above

can be evaluated by means of

which is the same result as before.

Gaussian quadrature#

A few functions are also provided in order to perform simple Gaussian quadrature over a fixed interval. The first is fixed_quad , which performs fixed-order Gaussian quadrature. The second function is quadrature , which performs Gaussian quadrature of multiple orders until the difference in the integral estimate is beneath some tolerance supplied by the user. These functions both use the module scipy.special.orthogonal , which can calculate the roots and quadrature weights of a large variety of orthogonal polynomials (the polynomials themselves are available as special functions returning instances of the polynomial class — e.g., special.legendre ).

Romberg Integration#

Romberg’s method [WPR] is another method for numerically evaluating an integral. See the help function for romberg for further details.

Integrating using Samples#

If the samples are equally-spaced and the number of samples available is \(2^+1\) for some integer \(k\) , then Romberg romb integration can be used to obtain high-precision estimates of the integral using the available samples. Romberg integration uses the trapezoid rule at step-sizes related by a power of two and then performs Richardson extrapolation on these estimates to approximate the integral with a higher degree of accuracy.

In case of arbitrary spaced samples, the two functions trapezoid and simpson are available. They are using Newton-Coates formulas of order 1 and 2 respectively to perform integration. The trapezoidal rule approximates the function as a straight line between adjacent points, while Simpson’s rule approximates the function between three adjacent points as a parabola.

For an odd number of samples that are equally spaced Simpson’s rule is exact if the function is a polynomial of order 3 or less. If the samples are not equally spaced, then the result is exact only if the function is a polynomial of order 2 or less.

How to solve differential equation using Python builtin function odeint?

I want to solve this differential equations with the given initial conditions:

the ans should be

here is my code:

but what I get is different from the answer. what have I done wrong?

1 Answer 1

There are several things wrong here. Firstly, your equation is apparently

(note the sign of the term in y). For this equation, your analytical solution and definition of y2 are correct.

Secondly, as the @Warren Weckesser says, you must pass 2 parameters as y to g : y[0] (y), y[1] (y’) and return their derivatives, y’ and y».

Thirdly, your initial conditions are given for x=0, but your x-grid to integrate on starts at -2. From the docs for odeint , this parameter, t in their call signature description:

odeint(func, y0, t, args=(). ) :

t : array A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.

So you must integrate starting at 0 or provide initial conditions starting at -2.

Finally, your range of integration covers a singularity at x=1/3. odeint may have a bad time here (but apparently doesn’t).