BTexample.mws

>    restart: assume(n,integer); assume(sigma,real); additionally(sigma>0);

>    bfn:=(x,n,a,b)->sin(n*Pi*(x+b)/(a+b))*sqrt(2/(a+b));

bfn := proc (x, n, a, b) options operator, arrow; sin(n*Pi*(x+b)/(a+b))*sqrt(2/(a+b)) end proc

>    BTLterm := (x,t,n,mu,sigma,a,b)->exp(mu*x/sigma^2-t*mu^2/(2*sigma^2))*bfn(0,n,a,b)*bfn(x,n,a,b)*exp((-t/2)*(n*Pi*sigma/(a+b))^2):

>    'BTLterm' = BTLterm(x,t,n,mu,sigma,a,b);

BTLterm = 2*exp(mu*x/sigma^2-1/2*t*mu^2/sigma^2)*sin(n*Pi*b/(a+b))/(a+b)*sin(n*Pi*(x+b)/(a+b))*exp(-1/2*t*n^2*Pi^2*sigma^2/(a+b)^2)

>    dBdt:=diff(BTLterm(x,t,n,mu,sigma,a,b),t): dBdxx:=diff(diff(BTLterm(x,t,n,mu,sigma,a,b),x),x): dBdx:=diff(BTLterm(x,t,n,mu,sigma,a,b),x):

>    ee:=dBdt - ((sigma^2/2)*dBdxx - mu*dBdx): 'PDE'=evalb(0=simplify(ee));

PDE = true

>   

>   

>    IBTLterm:=(int(BTLterm(x,1,n,mu,sigma,a,b),x=-b..a));

IBTLterm := -2*n*Pi*sigma^4*sin(n*Pi*b/(a+b))*exp(-1/2*(n^2*Pi^2*sigma^4+2*mu*b*a^2+4*mu*b^2*a+2*mu*b^3+mu^2*a^2+2*mu^2*a*b+mu^2*b^2)/sigma^2/(a+b)^2)*(-1+(-1)^n*exp(mu*(a+b)/sigma^2))/(mu^2*a^2+2*mu^2...

>    fhlcterm:=simplify(diff(diff(BTLterm(x,1,n,mu,sigma,a,b),a),b)):

>    fhlterm:=simplify(diff(diff(IBTLterm,a),b)):

>   

>    'fhlc(h,l,x,mu,sigma)'=(Sum(Diff(Diff(BTLterm(x,1,n,mu,sigma,a,b),a),b),n=1..infinity));

fhlc(h,l,x,mu,sigma) = Sum(Diff(2*exp(mu*x/sigma^2-1/2*mu^2/sigma^2)*sin(n*Pi*b/(a+b))/(a+b)*sin(n*Pi*(x+b)/(a+b))*exp(-1/2*n^2*Pi^2*sigma^2/(a+b)^2),a,b),n = 1 .. infinity)

>    'fhl(h,l,x,mu,sigma)'=Int(Sum(Diff(Diff(BTLterm(x,t,n,mu,sigma,a,b),a),b),n=1..infinity),x=-b..a);

fhl(h,l,x,mu,sigma) = Int(Sum(Diff(2*exp(mu*x/sigma^2-1/2*t*mu^2/sigma^2)*sin(n*Pi*b/(a+b))/(a+b)*sin(n*Pi*(x+b)/(a+b))*exp(-1/2*t*n^2*Pi^2*sigma^2/(a+b)^2),a,b),n = 1 .. infinity),x = -b .. a)

>    (Sum(Diff(Diff(IBTLterm,a),b),n=1..infinity));

Sum(Diff(-2*n*Pi*sigma^4*sin(n*Pi*b/(a+b))*exp(-1/2*(n^2*Pi^2*sigma^4+2*mu*b*a^2+4*mu*b^2*a+2*mu*b^3+mu^2*a^2+2*mu^2*a*b+mu^2*b^2)/sigma^2/(a+b)^2)*(-1+(-1)^n*exp(mu*(a+b)/sigma^2))/(mu^2*a^2+2*mu^2*a*...
Sum(Diff(-2*n*Pi*sigma^4*sin(n*Pi*b/(a+b))*exp(-1/2*(n^2*Pi^2*sigma^4+2*mu*b*a^2+4*mu*b^2*a+2*mu*b^3+mu^2*a^2+2*mu^2*a*b+mu^2*b^2)/sigma^2/(a+b)^2)*(-1+(-1)^n*exp(mu*(a+b)/sigma^2))/(mu^2*a^2+2*mu^2*a*...

--------------------------------------------------------------------------------

>    N:=16;

>    fhlc:=sum(fhlcterm,n=1..N): fhl:=sum(fhlterm,n=1..N):

>    ff:=(fhlc/fhl):

N := 16

>    aval:=1/2: bval:=1: muval:=0.0: sigval:=1:

>   

>    Gx:=evalf(subs(a=aval,b=bval,sigma=sigval,mu=muval,ff)):

>    int(Gx,x=-bval..aval);

.9999999995

>    plot(Gx,x=-bval..aval);

[Maple Plot]

>   

>   

>