[Python] numerical integration using smypy doesn't evaluate the expression

Discussão em 'Python' iniciado por Stack, Dezembro 22, 2024.

  1. Stack

    Stack Membro Participativo

    The integral is [​IMG]

    Where the variables follow the following distribution:

    [​IMG]

    Thus, the integral becomes:
    [​IMG]

    Now, my code in python is to integrate the expression with \(W_{1:3}\) being \(1,2,3\) respectively, \(r_1 = 0\), \(v\), which is the variance is \(1\), and the time constant \(tau\) being \(1\), I want to estimate roots of the derivative to find \(D\) as following:


    import sympy as sp
    import math

    w_n = [0,1,2]
    r_1 = sp.Symbol("r_1")
    r_2 = sp.Symbol("r_2")
    r_3 = sp.Symbol("r_3")
    r_n = [0, r_2,r_3]
    D = sp.Symbol("D")
    v = sp.Symbol("v")
    t = sp.Symbol("t")

    def expo_power(i=0,t=t,D=D,v=v):
    return ((((w_n - r_n) ** 2)/(2 * v)) + (((r_n - r_n[i-1]) ** 2)/(4 * D * t)))

    def expo_power_sum(lower_bound = 2, upper_bound=3,t=t,D=D,v=v):
    f = 0
    for i in range(lower_bound,upper_bound):
    f += expo_power(i-1,t=t,D=D,v=v)
    return f

    def P_w_n_P_r_n(i=0,v=v,D=D,t=t):
    return (1/(8 * (sp.pi ** 2) * D * t * v)) \
    * sp.exp(-expo_power_sum(lower_bound=2,upper_bound=3,t=t,D=D,v=v))

    def P_w_i_r_i(i = 0, v=v):
    return (1/(sp.sqrt(2 * sp.pi * v))) \
    * sp.exp(-((w_n - r_n) ** 2/(2 * v)))

    def normal_dist(x =0, m = r_n[0], v =v):
    return (1/(sp.sqrt(2 * sp.pi * v))) \
    * sp.exp(-(((x - m) ** 2)/(2 * v)))

    def integrand(v = v,t=t,lower = 2):
    f = P_w_n_P_r_n(i=lower,v=v,D=D,t=t)
    sp.pprint(f)
    f = f * P_w_i_r_i(v=v)
    f = f * normal_dist(x=r_n[0],v=v)
    return f

    def integrate(v=v,t=t, lower_bound = -sp.oo, upper_bound= sp.oo):
    function = integrand(v=v,t=t)
    integral_1 = sp.integrate(function, (r_1, lower_bound, upper_bound)).evalf()
    sp.pprint(integral_1)
    integral_2 = sp.integrate(integral_1, (r_3, lower_bound, upper_bound)).evalf()
    sp.pprint(integral_2)
    integral_3 = sp.integrate(integral_2, (r_2, lower_bound, upper_bound)).evalf()
    num = sp.N(integral_3)
    print(num)
    #sp.pprint(integral_3)
    deriv = sp.diff(num, D)
    sol = sp.solve(deriv, D)
    sp.pprint(sol)
    integrate(v=1,t=1,lower_bound=-1000,upper_bound=1000)



    Now, the solution displayed is the ratio between two unevaluated integrals.

    note that [​IMG] and t in the python code is (tau) and we're finding the roots of the derivative with respect to (D) after marginalizing out (r_1,r_2,r_3), the integral is not being evaluated, particularly, for (r_2)

    Continue reading...

Compartilhe esta Página