How to solve the derivative of the expression with respect to D?

The integral is

Where the variables follow the following distribution:

Thus, the integral becomes:

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

N = 3

w_n = [k for k in range(N)]
r_n = [sp.symbols(f'r_{i+1}') for i in range(N)]
D = sp.symbols("D",real=True,)
v = sp.symbols("v",real=True)
t = sp.symbols("t",real=True)

def expo_power(i=0,t=t,D=D,v=v):
    f = ((((w_n[i] - r_n[i]) ** 2)/(2 * v)) + (((r_n[i] - r_n[i-1]) ** 2)/(4 * D * t)))
    return f
def expo_power_sum(lower_bound = 2, upper_bound=N,t=t,D=D,v=v):
    f = 0
    for i in range(lower_bound, upper_bound):
        f += expo_power(i,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/sp.sqrt(8 * (sp.pi ** 2) * D * t * v) ** (N - 1) ) \
        * sp.exp(-expo_power_sum(lower_bound=2,upper_bound=N,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[i] - r_n[i]) ** 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 = sp.log(P_w_n_P_r_n(i=lower,v=v,D=D,t=t) * P_w_i_r_i(v=v) * 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)

    for i in range(N):  
        inte = sp.Integral(function,(r_n[i],lower_bound,upper_bound))
        inte = inte.doit()
        inte = inte.evalf()
        function = inte
    k = inte
    sp.pprint(k)

    d = sp.diff(k,D)
    sol = sp.solve(d,D)
    sp.pprint(sol)
integrate(v=1,t=1,lower_bound=0,upper_bound=10)

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

note that 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)

Does the integral have a closed form? if so can the derivative of the form be solved for $D$? if not, what are the alternative to produce a solution for $P'(w_{1:3}|D,v) = 0$ with respect to $D$?

EDIT: code edited after fixing bugs pointed out by @RoberDodier

I can't tell fully what's going on, but there seem to be some problems with the way the problem is stated, and I'm thinking that you need to straighten out the problem statement before seeking a solution.

If I'm not mistaken, all the distributions involved are normal distributions, and you are marginalizing (integrating) over some of them. If so, the result should again be a normal distribution (if there are still some free variables) or some constant involving the mean and variance. My guess is that if you phrased the joint distribution as a multivariate normal distribution with a suitable covariance matrix, you would be able to state the result in terms of operations on the covariance.

But even without that, it should be possible to get an exact result in terms of D, v, and t, if the problem is set up right.

The first problem I see is that the integrand (function in your function integrate) is a function of r_2 only. It should be a function of r_1 and r_3 as well, shouldn't it?

Shouldn't r_n be assigned [r_1, r_2, r_3] and not [0, r_2, r_3] ?

That function is free of r_3 suggests that expo_power_sum isn't working right.

Finally, my advice is to work with symbolic infinity (plus/minussp.oo) when stating the integrals instead of plus/minus 1000 -- substituting any finite value is likely to make the results more complicated than they need to be.