Как решить интегральное уравнение в питоне? - PullRequest
Новая
       40

Как решить интегральное уравнение в питоне?

3 голосов
/ 23 августа 2011

Я пытаюсь решить это интегральное уравнение, используя Python:

enter image description here

, где z колеблется от 0 до 1.

Функция scipy.quad предоставляет численное решение только для определенного интервала, но не обеспечивает решение в течение интервала. 

def f(z,Om,Ol): return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
quad(lambda r:f(r,Om,Ol),0,1)
(0.77142706642781111, 8.5645609096719596e-15)

Но я не знаю, как получить полный вектор в этом интервале, как вы получаете, когда используете scipy.odeint для решения дифференциального уравнения.

С другой стороны, sympy.integrate не может этого сделать. Я получаю переполнение стека. Кроме того, я не могу понять, как заменить символы списком, то есть: 

sy.integrate(x**2,x).subs(x,1)
1/3
sy.integrate(x**2,x).subs(x,[1,2])
TypeError: unhashable type: 'list'

Так что вопрос: кто-нибудь знает, как решить интегральное уравнение, используя python?

Ответы [ 4 ]

4 голосов
/ 23 августа 2011

Я понимаю, что вы хотите решить дифференциальное уравнение dF/dz1 = f(z1, Om, Ol) и хотите F(z1) в разных местах.Если это так, то подпрограмма Обыкновенного дифференциального уравнения (ODE) *1004* SciPy - это то, что нужно.В частности, вы можете проверить odeint(), так как он может дать вам значения вашего интеграла в указанных вами местах.

1 голос
/ 23 августа 2011

Я полагаю, что z перед интегралом является опечаткой, которая должна быть z1, и вы ищете z1 с учетом DL.

Сначала вы должны реализовать правую часть (rhs): 

def f(z,Om,Ol): 
        return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)
def rhs(z1, Om, Ol, c, H0):
        return c/H0*(1+z1)*quad(lambda r:f(r, Om, Ol), 0, z1)[0]

Теперь вам нужно найти z0 такое, что rhs (z1, ...) = DL, что совпадает с 

rhs(z1, ...) - DL = 0

Это означает, что ваша проблема сводится к нахождению нуля (есть только один, потому что rhs является монотонным), 

f(z1) = rhs(z1, ...) - DL

Здесь вы можете применить множество методов для поиска корней (см. http://en.wikipedia.org/wiki/Root-finding_algorithm) из http://docs.scipy.org/doc/scipy/reference/optimize.html#root-finding

0 голосов
/ 26 августа 2011

Я наконец-то использовал функцию "quad" с оператором for для решения моей проблемы: 

import pylab as p
import scipy as s
from scipy.integrate import odeint,quad

def Dl_lflrw(z,Om,Ol):
    c = 1.
    H0 = 1.
    y = []
    for i in z:
        y1 = (c/H0)*(1+i)*quad(lambda r:f(r,Om,Ol),0,i)[0]
        y.append(y1)
    return p.asarray(y)

def f(z,Om,Ol):
    return 1./p.sqrt((1+z)**2 * (1+Om*z) - z*(2+z)*Ol)

Спасибо всем за ваши идеи, они действительно помогли.

0 голосов
/ 23 августа 2011

В разделе примеров в разделе http://docs.sympy.org/0.7.1/modules/integrals.html, показаны решения практически идентичных проблем.Можете ли вы опубликовать свой код Simpy?

Для scipy, вы пытались использовать кортеж, который можно хэшировать, вместо списка?например: 

sy.integrate(x**2,x).subs(x,(1,2,))
...