We are going to show a roundabout but rather effective method of solving some rather complicated (definite) integrals in closed-form via a rather unusual method: a special factorization of linear ODEs.  The example we will use is a 2 week old question that has yet to get an answer.

First, the problem: compute the integral

f  := exp(a*m)* (sqrt(T)*sqrt(Pi)*erf((1/2)*sqrt(2)*(m+a*T)/sqrt(T))*exp(a*m)*a -sqrt(T)*sqrt(Pi)*exp(a*m)*a +sqrt(2)*exp(-(1/2)*(m^2+a^2*T^2)/T)):  gg := Int(f,m=0..K);

Int(exp(a*m)*(T^(1/2)*Pi^(1/2)*erf(1/2*2^(1/2)*(m+a*T)/T^(1/2))*exp(a*m)*a -T^(1/2)*Pi^(1/2)*exp(a*m)*a+2^(1/2)*exp(-1/2*(m^2+a^2*T^2)/T)),m = 0 .. K)

The first step is to translate this into a problem on ODEs:

dd := PDEtools[dpolyform](y(K)=gg,no_Fn);
dd2 := numer((lhs-rhs)(op([1,1],dd)));

which gives a 4th order linear ODE with polynomial coefficients. Then we need to draw on the power of the 'diffop' (for differential operators) sub-package of DEtools, in particular the LCLM factorization.

oper := DEtools[de2diffop](dd2, y(K), [Dx,x]);
fact := DEtools[DFactorLCLM](oper,[Dx,x]);
de3 := map(DEtools[diffop2de], fact, [Dx,x], y(x));

[-(-x+a*T)/T*diff(y(x),x)+diff(diff(y(x),x),x), 2*(-x+a*T)/T*a*y(x)-1/T*(-x+3*a*T)*diff(y(x),x)+diff(diff(y(x),x),x)]

So now we have something quite fantastic: the solution to the integral will be made up of the solutions to those 2 second-order ODEs. What will those look like?

basis4  := map(op,map(dsolve,de3,output=basis));
eq := add(C[i]*basis4[i],i=1..4);

C[1]+C[2]*erf((-1/2*x+1/2*a*T)/T^(1/2)*2^(1/2))+ C[3]*exp(2*a*x)+C[4]*exp(2*a*x)*erf((1/2*a*T+1/2*x)*2^(1/2)/T^(1/2))

Those are exactly as asked by the original poster! We need to find the exact solution, so we just form the linear system from the first 4 Taylor coefficients, solve it, and substitute

eqs := [seq(eval(diff(eq,[x$i]),x=0) - eval(value(diff(gg,[K$i])),K=0),i=0..3)]:
solve(eqs,{C[1],C[2],C[3],C[4]});
eval(eq,%);

1/2*Pi^(1/2)*T^(1/2)-1/2*Pi^(1/2)*T^(1/2)*erf((-1/2*x+1/2*a*T)/T^(1/2)*2^(1/2))-1/2*Pi^(1/2)*T^(1/2)*exp(2*a*x)+1/2*Pi^(1/2)*T^(1/2)*exp(2*a*x)*erf((1/2*a*T+1/2*x)*2^(1/2)/T^(1/2))


Please Wait...