#COPYRIGHT NOTICE: Copyright (c) 1995-98 by Sasha Cyganowski. Permission is
#granted to anyone to use, modify, or redistribute this software freely,
#subject to the following restrictions: 1. The author accepts no
#responsibility for any consequences of this software and makes no guarantee
#that this software is free of defects. 2. The origin of this software must
#not be misrepresented, either by explicit claim or by omission.
#3. This notice and the copyright must be included with all copies or
#modified versions of this software. 4. This software may not be included or
#redistributed as part of any package to be sold for profit unlesss the author
#has given explicit written permission to do so.
# Sasha Cyganowski, School of Computing and Mathematics, Deakin University,
#Waurn Ponds, VIC, Australia. Internet: sash@deakin.edu.au
# Telephone: +61 3 5227 1382.
# This file contains the function definitions and help files for the package
#"stochastic", developed by Sasha Cyganowski. The package contains functions
#which can be used to find explicit solutions of Stochastic Differential
#Equations (SDEs), construct numerical schemes up to strong order 2 and weak
#order 3, check for commutative noise of the first and second kind, and
#convert SDEs into their coloured noise form. Other useful routines are also
#included.
#The general schemes, solutions, and conditions from which the code
#was developed can be found in Kloeden, P.E, Platen, E.: Numerical
#Solution of Stochastic Differential Equations, (Springer-Verlag, 1992).
#
#
#A manual for this package is available from
#
# http://www.math.uni-frankfurt.de/~numerik/maplestoch/
#
#This version of the software is known to work for Maple V, Release 5.1.
#
#Additional Routines by Thomas Pohl
#
#Updates and bugfixes by Lars Gruene
#http://www.uni-bayreuth.de/departments/math/~lgruene/
#
#Last Change: December 11, 2003 by Lars Gruene
`help/text/stochastic` := TEXT(
" ",
"HELP FOR: Introduction to the stochastic package",
" ",
"CALLING SEQUENCE:",
" (args)",
" stochastic[](args)",
" ",
"SYNOPSIS:",
"- The stochastic package contains routines useful for finding explicit",
" solutions of Stochastic Differential Equations (SDEs), and routines",
" useful for constructing numerical schemes up to strong order 2 and",
" weak order 3.",
" Other useful features include a routine which converts SDEs with white",
" noise into coloured noise form, and routines which check whether",
" an SDE has commutative noise of the first and/or second kind. Other",
" useful routines are also available. In particular, the operator",
" procedures L0 and LJ are included so that users can easily construct",
" numerical schemes other than those already available.",
" For more information on a particular function type the command",
" ?.",
" The author of the stochastic package is Sasha Cyganowski, from Deakin",
" University, Waurn Ponds, Australia, 3217. Comments are welcomed via",
" e-mail, which can be sent to the folllowing address: sash@deakin.edu.au",
" The general schemes, solutions, and conditions used in the coding of the",
" stochastic package can be found in Kloeden, P.E, Platen, E.: Numerical",
" Solution of Stochastic Differential Equations, (Springer-Verlag, 1992).",
" The package is also described in Cyganowski, S.O., A MAPLE Package for",
" Stochastic Differential Equations, Computational Techniques and",
" Applications: CTAC95, editors A. Easton, R. May, World Scientific,",
" (to appear).",
" ",
"- To use a stochastic function, either define that function alone using the",
" command with(stochastic, ), or define all stochastic functions",
" using the command with(stochastic). Alternatively, invoke the function",
" using the long form stochastic[]. This long form notation is",
" necessary whenever there is a conflict between a package function name",
" and another function used in the same session.",
" ",
"- The functions available are:",
" ",
" linearsde MLJ explicit L0 Euler Milstein",
" Taylor1hlf SLO correct Taylor2 wkeuler wktay2",
" reducible wktay3 colour comm1 comm2 LJ",
" ",
"- As an example, to find the explicit solution of an SDE with drift",
" coefficient 1/2*a^2*x and diffusion coefficient a*x, use",
" ",
" with(stochastic,explicit); explicit(1/2*a^2*x,a*x);",
" ",
"- Note the stochastic numerical routines require the drift and diffusion",
" coefficients of an SDE to be entered in the variables x[N] and t, where",
" x[N] are the state variables of the N-dimensional SDE and t denotes time.",
" The routines linearsde, reducible, and explicit require the drift and",
" duffusion coefficients to be entered in the variables x and t.",
" ",
"SEE ALSO: with"
):
stochastic[linearsde]:=proc(a::algebraic,b::algebraic)
local temp1,alpha,beta,gamma,delta,fundsoln,fundsoln2,soln,default1,default2,
default3;
if diff(a,x,x) <> 0 or diff(b,x,x) <> 0 then
ERROR(`SDE not linearsde, try a reducible procedure`)
else
alpha := diff(a,x);
alpha := subs(t = s,alpha);
beta := diff(b,x);
beta := subs(t = s,beta);
if diff(beta,s) = 0 then temp1 := beta*W;
else temp1:=Int(beta,W = 0 .. t);
fi;
gamma := coeff(a,x,0);
gamma := subs(t = s,gamma);
delta := coeff(b,x,0);
delta := subs(t = s,delta);
fundsoln := exp(int(alpha-1/2*beta^2,s = 0 .. t)+temp1);
fundsoln2 := subs(t = s,fundsoln);
if beta = 0 then
soln := fundsoln*(X[0]+int(1/fundsoln2*(gamma-beta*delta),
s = 0 .. t)+Int(1/fundsoln2*delta,W = 0 .. t))
else soln := fundsoln*(X[0]+Int(1/fundsoln2*(gamma-beta*delta),
s = 0 .. t)+Int(1/fundsoln2*delta,W = 0 .. t))
fi;
default1 := Int(0,W = 0 .. t) = 0;
default2 := Int(0,W = 0 .. s) = 0;
default3 := Int(0,s = 0 .. t) = 0;
soln := X[t] = subs(default1,default2,default3,soln);
fi;
end:
`help/text/linearsde` := TEXT(
"FUNCTION: stochastic[linearsde] - returns the explicit solution of a linearsde",
"Stochastic Differential Equation (SDE).",
" ",
"CALLING SEQUENCE: linearsde(a,b);",
" ",
"PARAMETERS: a - algebraic, given in the variables x and t.",
" b - algebraic, given in the variables x and t.",
" ",
"SYNOPSIS: ",
"- The call linearsde(a,b); returns the explicit solution for a linearsde",
" Stochastic Differential Equation (SDE) with drift coefficient a, and",
" diffusion coefficient b. The output consists of the variables X[t], X[0],",
" and W. X[t] denotes the explicit solution, X[0] denotes the initial value",
" of the solution, W denotes a standard Wiener process, and t denotes time.",
" If the SDE is not of linearsde form a suitable error message is returned.",
" ",
"- This routine is used by the routine stochastic[explicit] and, in general,",
" is not intended to be used on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[linearsde].",
" ",
"EXAMPLES: ",
" ",
"linearsde(-x,2);",
" ",
" / t \ ",
" | / |",
" | | 2 |",
" X[t] = exp(- t) |X[0] + | -------- dW|",
" | | exp(- s) |",
" | / |",
" \ 0 /",
" ",
"SEE ALSO: stochastic, stochastic[reducible], stochastic[explicit] "
):
stochastic[reducible]:=proc(a::algebraic,b::algebraic)
local beta,temp1,h,temp3,alpha,soln,soln1,II;
h := int(1/b,x);
temp1 := alpha*b*h+1/2*b*simplify(diff(b,x));
temp1 = a;
alpha := simplify(solve(%,alpha));
beta := alpha*h;
if diff(alpha,x) = 0 then
if alpha=0 then
soln:=h=subs(x=X[0],h)+W;
X[t]=simplify(solve(soln,x));
else
soln1 := h = exp(alpha*t)*subs(x = X[0],h)+exp(alpha*t)*II;
X[t] = subs(II=Int(exp(-alpha*s),W=0..t),solve(soln1,x));
fi
elif diff(beta,x) = 0 then
X[t]=simplify(solve(h = beta*t+W+subs(x = X[0],h),x));
else ERROR(`non-linearsde SDE not reducible`)
fi
end:
`help/text/reducible` := TEXT(
"FUNCTION: stochastic[reducible] - returns the explicit solution of a",
"reducible Stochastic Differential Equation (SDE).",
" ",
"CALLING SEQUENCE: reducible(a,b);",
" ",
"PARAMETERS: a - algebraic, given in the variables x and t.",
" b - algebraic, given in the variables x and t.",
" ",
"SYNOPSIS: ",
"- The call reducible(a,b); returns the explicit solution for a reducible",
" Stochastic Differential Equation (SDE) with drift coefficient a, and",
" diffusion coefficient b. The output consists of the variables X[t], X[0],",
" and W. X[t] denotes the explicit solution, X[0] denotes the initial",
" value of the solution, W denotes a standard Wiener process, and t denotes",
" time. If the SDE is not of reducible form a suitable error message is",
" returned.",
" ",
"- This routine is used by the routine stochastic[explicit] and, in general,",
" is not intended to be used on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be called via the call",
" stochastic[reducible].",
" ",
"EXAMPLES:",
" ",
">reducible(1,2*sqrt(x));",
" ",
" 1/2 2",
" X[t] = 2 X[0] W + W + X[0]",
" ",
"SEE ALSO: stochastic, stochastic[linearsde],stochastic[explicit] "
):
stochastic[explicit]:=
proc(a::algebraic,b::algebraic) if diff(a,x,x) = 0 and diff(b,x,x) = 0 then linearsde
(a,b) else reducible(a,b) fi end:
`help/text/explicit` := TEXT(
"FUNCTION: stochastic[explicit] - returns the explicit solution of a",
" Stochastic Differential Equation (SDE).",
" ",
"CALLING SEQUENCE: explicit(a,b);",
" ",
"PARAMETERS: a - algebraic, given in the variables x and t.",
" b - algebraic, given in the variables x and t.",
" ",
"SYNOPSIS: ",
"- The call explicit(a,b); returns the explicit solution for a Stochastic",
" Differential Equation (SDE) with drift coefficient a, and diffusion",
" coefficient b. The output consists of the variables X[t], X[0], and W.",
" X[t] denotes the explicit solution, X[0] denotes the initial value of the",
" solution, and W denotes a standard Wiener process. The user is returned",
" with a suitable error message if no known explicit solution exists.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[explicit].",
" ",
"EXAMPLES: ",
" ",
">explicit(1/2*a^2*x,a*x);",
" ",
" X[t]=exp(aW)X[0]",
" ",
"SEE ALSO: stochastic, stochastis[linearsde], stochastic[reducible]"
):
stochastic[LJ]:=
proc(X::algebraic,b::list(list(algebraic)),j::integer) sum('op(j,op(k,b))*
diff(X,x[k])','k' = 1 .. nops(b)) end:
`help/text/LJ` := TEXT(
"FUNCTION: stochastic[LJ] - applies the partial differential operator LJ, as",
" described in Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic",
" Differential Equations, pg 339, (Springer-Verlag 1992).",
" ",
"CALLING SEQUENCE: LJ(X,[[b11,..,b1M],..,[bN1,..,bNM]],j);",
" ",
"PARAMETERS: X - algebraic, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" j - integer.",
" ",
"SYNOPSIS: ",
"- The call LJ(X,[[b11,..,b1M],..,[bN1,..,bNM]],j); computes the application",
" of the partial differential operator LJ to X, as described in",
" Kloeden, P.E., Platen, E: Numerical Solution of Stochastic Differential",
" Equations, pg 339, (Springer-Verlag 1992).",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the",
" diffusion matrix of dimension N*M, where M is the dimension of the",
" Wiener process and N is the dimension of the SDE. j=1,..,M denotes the",
" ''current'' dimension of the Wiener process.",
" The output variables are consistent with the variables used as input.",
" ",
"- This routine is used by the following routines: stochastic[Euler],",
" stochastic[Milstein], stochastic[Taylor1.5], stochastic[Taylor2],",
" stochastic[wkeuler], and stochastic[wktay2].",
" In general, it is not intended for use on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[LJ]",
" ",
"EXAMPLES: ",
" ",
">LJ(x[2],[[r,0],[0,r]],2);",
" ",
" r ",
" ",
">LJ(x[2],[[r,0],[0,r]],1);",
" ",
" 0 ",
" ",
"SEE ALSO: stochastic[L0], stochastic[MLJ], stochastic[SL0],",
" stochastic[Euler], stochastic[Milstein], stochastic[Taylor1hlf]",
" stochastic[Taylor2], stochastic[wkeuler], stochastic[wktay2],",
" stochastic."
):
stochastic[L0]:=proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2,part3;
part1 := diff(X,t);
part2 := sum('a[k]*diff(X,x[k])','k' = 1 .. nops(a));
part3 := 1/2*sum(
'sum('sum('op(j,op(k,b))*op(j,op(l,b))*diff(X,x[k],x[l])',
'j' = 1 .. nops(op(1,b)))','k' = 1 .. nops(a))','l' = 1 .. nops(a)
);
part1+part2+part3;
end:
`help/text/L0` := TEXT(
"FUNCTION: stochastic[L0] - applies the partial differential operator L0, as",
" described in Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic",
" Differential Equations, pg 339, (Springer-Verlag 1992).",
" ",
"CALLING SEQUENCE: L0(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: X - algebraic, given in the variables x[N] and t.",
" a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call L0(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); computes the",
" application of the partial differential operator L0 to X, as described",
" in Kloeden, P.E., Platen, E: Numerical Solution of Stochastic",
" Differential Equations, pg 339, (Springer-Verlag 1992). a1,..,aN",
" denotes the drift coefficients of the N-dimensional SDE.",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M.",
" The output variables are consistent with the variables used as input.",
" ",
"- This routine is used by the following routines: stochastic[Euler],",
" stochastic[Milstein], stochastic[Taylor1.5], stochastic[wkeuler]",
" and stochastic[wktay2].",
" In general, it is not intended for use on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[L0]",
" ",
"EXAMPLES: ",
" ",
">L0(x[2],[x[1],x[2]],[[r,0],[0,r]]);",
" ",
" x[2] ",
" ",
"SEE ALSO: stochastic[LJ], stochastic[MLJ], stochastic[SL0],",
" stochastic[Euler], stochastic[Milstein], stochastic[Taylor1hlf]",
" stochastic[wkeuler], stochastic[wktay2], stochastic[wktay3],",
" stochastic."
):
stochastic[Euler]:=proc(a::list(algebraic),b::list(list(algebraic)))
local i,u,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+L0(x[i],a,b)*Delta[n]+sum('LJ(x[i],b,j)*
Delta*W.j[n]','j' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/Euler` := TEXT(
"FUNCTION: stochastic[Euler] - constructs stochastic Taylor schemes of",
"strong order 0.5, otherwise known as Euler schemes.",
" ",
"CALLING SEQUENCE: Euler([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call Euler([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" Euler scheme for an N-dimensional Stochastic Differential Equation (SDE)",
" with M-dimensional noise which has as its drift coefficients a1,..,aN and",
" [b11,..,b1M],..,[bN1,..,bNM] as its diffusion coefficients.",
" The output consists of the variables YN[n], DeltaWM[n], and Delta[n].",
" YN[n] denotes the strong order 0.5 stochastic Taylor approximation",
" to x[N] at the n-th step. DeltaWM[n] denotes the change in the",
" M-dimensional Wiener process at the n-th step. Delta[n] denotes the step",
" size at the n-th step.",
" ",
" This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[Euler].",
" ",
"EXAMPLES: ",
" ",
">Euler([x[2],x[1]],[[r,s],[s,r]]);",
" ",
"table([",
" 1 = (Y1[n + 1] = Y1[n] + Y2[n] Delta[n] + r Delta W1[n] + s Delta W2[n])",
" 2 = (Y2[n + 1] = Y2[n] + Y1[n] Delta[n] + s Delta W1[n] + r Delta W2[n])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ], ",
" stochastic[Milstein], stochastic[Taylor1hlf], stochastic[Taylor2]"
):
stochastic[Milstein]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+L0(x[i],a,b)*Delta[n]+sum('LJ(x[i],b,j)*
Delta*W.j[n]','j' = 1 .. nops(op(1,b)))+
sum('sum('LJ(op(j2,op(i,b)),b,j1)*I[j1,j2]',
'j1' = 1 .. nops(op(1,b)))','j2' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/Milstein` := TEXT(
"FUNCTION: stochastic[Milstein] - constructs stochastic Taylor schemes of",
"strong order 1, otherwise known as Milstein schemes.",
" ",
"CALLING SEQUENCE: Milstein([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call Milstein([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" Milstein scheme for an N-dimensional Stochastic Differential Equation",
" (SDE) with M-dimensional noise which has as its drift coefficients",
" a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its diffusion coefficients.",
" The output consists of the variables YN[n], DeltaWM[n], Delta[n],",
" and I[(j1,j2)]. ",
" YN[n] denotes the strong order 1 stochastic Taylor approximation",
" to x[N] at the n-th step. DeltaWM[n] denotes the change in the",
" M-dimensional Wiener process at the n-th step. Delta[n] denotes the step",
" size at the n-th step. I[(j1,j2)] denotes multiple Ito integrals",
" (refer to Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic",
" Differential Equations, (Springer-Verlag, 1992)).",
" ",
" This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[Milstein].",
" ",
"EXAMPLES: ",
" ",
">Milstein([x[2],x[1]],[[r,s],[s,r]]);",
" ",
"table([",
" 1 = (Y1[n + 1] = Y1[n] + Y2[n] Delta[n] + r Delta W1[n] + s Delta W2[n])",
" 2 = (Y2[n + 1] = Y2[n] + Y1[n] Delta[n] + s Delta W1[n] + r Delta W2[n])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ],",
" stochastic[Euler], stochastic[Taylor1hlf], stochastic[Taylor2]"
):
stochastic[Taylor1hlf]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+a[i]*Delta[n]+1/2*L0(a[i],a,b)*Delta[n]^2+
sum('op(j,op(i,b))*Delta*W.j[n]+L0(op(j,op(i,b)),a,b)*I[0,j]+
LJ(a[i],b,j)*I[j,0]','j' = 1 .. nops(op(1,b)))+
sum('sum('LJ(op(j2,op(i,b)),b,j1)*I[j1,j2]',
'j1' = 1 .. nops(op(1,b)))','j2' = 1 .. nops(op(1,b)))+sum(
'sum('sum('LJ(LJ(op(p3,op(i,b)),b,p2),b,p1)*I[p1,p2,p3]',
'p1' = 1 .. nops(op(1,b)))','p2' = 1 .. nops(op(1,b)))',
'p3' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/Taylor1hlf` := TEXT(
"FUNCTION: stochastic[Taylor1hlf] - constructs stochastic Taylor schemes of",
"strong order 1.5.",
" ",
"CALLING SEQUENCE: Taylor1hlf([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call Taylor1hlf([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns",
" the strong order 1.5 approximation for an N-dimensional Stochastic",
" Differential Equation (SDE) with M-dimensional noise which has as its",
" drift coefficients a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its",
" diffusion coefficients.",
" The output consists of the variables YN[n], DeltaWM[n], Delta[n],",
" I[(j1,j2)], and I[(j1,j2,j3)]. ",
" YN[n] denotes the strong order 1.5 stochastic Taylor approximation",
" to x[N] at the n-th step. DeltaWM[n] denotes the change in the",
" M-dimensional Wiener process at the n-th step. Delta[n] denotes the step",
" size at the n-th step. I[(j1,j2)] and I[(j1,j2,j3)] denote multiple Ito",
" integrals (refer to Kloeden, P.E., Platen, E.: Numerical Solution of",
" Stochastic Differential Equations, (Springer-Verlag, 1992)).",
" ",
" This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[Taylor1hlf].",
" ",
"EXAMPLES: ",
" ",
">Taylor1hlf([x[2],x[1]],[[r,s],[s,r]]);",
" ",
"table([",
" 2",
" 1 = (Y1[n + 1] = Y1[n] + Y2[n] Delta[n] + 1/2 Y1[n] Delta[n] ",
" + r Delta W1[n] ",
" + s I[1, 0] + s Delta W2[n] + r I[2, 0])",
" 2",
" 2 = (Y2[n + 1] = Y2[n] + Y1[n] Delta[n] + 1/2 Y2[n] Delta[n] ",
" + s Delta W1[n] ",
" + r I[1, 0] + r Delta W2[n] + s I[2, 0])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ],",
" stochastic[Euler], stochastic[Milstein],",
" stochastic[Taylor2]"
):
stochastic[SL0]:=proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2;
part1 := diff(X,t); part2 := sum('a[k]*diff(X,x[k])','k' = 1 .. nops(a));
part1+part2;
end:
`help/text/SL0` := TEXT(
"FUNCTION: stochastic[SL0] - applies the partial differential operator (1.2),",
" as described in Kloeden, P.E., Platen, E: Numerical Solution of Stochastic",
" Differential Equations, pg 339, (Springer-Verlag 1992). This is the",
" Stratonovich version of the operator L0.",
" ",
"CALLING SEQUENCE: SL0(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: X - algebraic, given in the variables x[N] and t.",
" a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call SL0(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); computes the",
" application of the Stratonovich version of the operator L0 to X, as",
" described in Kloeden, P.E., Platen, E: Numerical Solution of Stochastic",
" Differential Equations, pg 339, (Springer-Verlag 1992).",
" a1,..,aN denotes the drift coefficients of the N-dimensional SDE.",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M, where M is the dimension of the Wiener process.",
" The output variables are consistent with the variables used as input.",
" ",
"- This routine is used by the routine stochastic[Taylor2].",
" In general, it is not intended for use on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[SL0]",
" ",
"EXAMPLES: ",
" ",
">SL0(x[2],[x[1],x[2]],[[r,0],[0,r]]);",
" ",
" x[2] ",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ],",
" stochastic[MLJ], stochastic[Taylor2]"
):
stochastic[correct]:=proc(a::list(algebraic),b::list(list(algebraic)),i) a[i]-
1/2*sum('LJ(op(j,op(i,b)),b,j)','j' = 1 .. nops(op(1,b))); end:
`help/text/correct` := TEXT(
"FUNCTION: stochastic[correct] - applies the drift-correction formula for",
"converting Ito Stochastic Differential Equations (SDEs) to Stratonovich",
"SDEs.",
" ",
"CALLING SEQUENCE: correct([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]],i);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" i - integer.",
" ",
"SYNOPSIS: ",
"- The call correct([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]],i); converts",
" the drift coefficient of an Ito SDE to its Stratonovich form.",
" The drift-correction formula can be found in Kloeden, P.E.,",
" Platen, E.: Numerical Solution of Stochastic Differential Equations,",
" pg 339, (Springer-Verlag, 1992).",
" a1,..,aN denotes the drift coefficients of the N-dimensional SDE.",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M, where M is the dimension of the Wiener process. i = 1..N denotes",
" the ''current'' dimension of the SDE.",
" The output variables are consistent with the variables used as input.",
" ",
"- This routine is used by the routine stochastic[Taylor2].",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[correct].",
" ",
"EXAMPLES: ",
" ",
">correct([x[1],x[2]],[[r,0],[0,r]],2);",
" ",
" x[2] ",
" ",
">correct([x[1],x[2]],[[x[1],0],[0,r]],1);",
" ",
" 1/2 x[1]",
" ",
"SEE ALSO: stochastic, stochastic[Taylor2]"
):
stochastic[Taylor2]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+correct(a,b,i)*Delta[n]+
1/2*SL0(correct(a,b,i),a,b)*Delta[n]^2+
sum('op(j,op(i,b))*Delta*W.j[n]+SL0(op(j,op(i,b)),a,b)*J[0,j]+
LJ(correct(a,b,i),b,j)*J[j,0]','j' = 1 .. nops(op(1,b)))
+sum('sum('LJ(op(j2,op(i,b)),b,j1)*J[j1,j2]+
SL0(LJ(op(j2,op(i,b)),b,j1),a,b)*J[0,j1,j2]+
LJ(SL0(op(j2,op(i,b)),a,b),b,j1)*J[j1,0,j2]+
LJ(LJ(correct(a,b,i),b,j2),b,j1)*J[j1,j2,0]',
'j1' = 1 .. nops(op(1,b)))',
'j2' = 1 .. nops(op(1,b)))+sum(
'sum('sum('LJ(LJ(op(p3,op(i,b)),b,p2),b,p1)*J[p1,p2,p3]',
'p1' = 1 .. nops(op(1,b)))','p2' = 1 .. nops(op(1,b)))',
'p3' = 1 .. nops(op(1,b)))+sum('sum('sum(
'sum('LJ(LJ(LJ(op(m4,op(i,b)),b,m3),b,m2),b,m1)*J[m1,m2,m3,m4]',
'm1' = 1 .. nops(op(1,b)))','m2' = 1 .. nops(op(1,b)))
','m3' = 1 .. nops(op(1,b)))','m4' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/Taylor2` := TEXT(
"FUNCTION: stochastic[Taylor2] - constructs stochastic Taylor schemes of",
"strong order 2.",
" ",
"CALLING SEQUENCE: Taylor2([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call Taylor2([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" strong order 2 stochastic Taylor approximation for an N-dimensional",
" Stochastic Differential Equation (SDE) with M-dimensional noise which has",
" as its drift coefficients a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its",
" diffusion coefficients.",
" The output consists of the variables YN[n], DeltaWM[n], Delta[n],",
" J[(j1,j2)], J[(j1,j2,j3)], and J[(j1,j2,j3,j4)].",
" YN[n] denotes the strong order 2 stochastic Taylor approximation",
" to x[N] at the n-th step. DeltaWM[n] denotes the change in the",
" M-dimensional Wiener process at the n-th step. Delta[n] denotes the step",
" size at the n-th step. J[(j1,j2)], J[(j1,j2,j3)], and J[(j1,j2,j3,j4)]",
" denote multiple Stratonovich integrals (refer Kloeden, P.E., Platen, E.:",
" Numerical Solution of Stochastic Differential Equations,",
" (Springer-Verlag, 1992). ",
" ",
" This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[Taylor2].",
" ",
"EXAMPLES: ",
" ",
">Taylor2([x[2],x[1]],[[r,s],[s,r]]);",
" ",
"table([",
" 2",
" 1 = (Y1[n + 1] = Y1[n] + Y2[n] Delta[n] + 1/2 Y1[n] Delta[n] ",
" + r Delta W1[n] ",
" + s J[1, 0] + s Delta W2[n] + r J[2, 0])",
" 2",
" 2 = (Y2[n + 1] = Y2[n] + Y1[n] Delta[n] + 1/2 Y2[n] Delta[n] ",
" + s Delta W1[n] ",
" + r J[1, 0] + r Delta W2[n] + s J[2, 0])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[LJ], stochastic[SL0],",
" stochastic[Euler], stochastic[Milstein], stochastic[Taylor1hlf],",
" stochastic[correct]."
):
stochastic[wkeuler]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+L0(x[i],a,b)*Delta[n]+
sum('LJ(x[i],b,j)*Delta*Ws.j[n]','j' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/wkeuler` := TEXT(
"FUNCTION: stochastic[wkeuler] - constructs stochastic Taylor schemes of",
"weak order 1, otherwise known as weak Euler schemes.",
" ",
"CALLING SEQUENCE: wkeuler([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call wkeuler([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" simplified weak euler scheme for an N-dimensional Stochastic Differential",
" Equation (SDE) with M-dimensional noise which has as its drift",
" coefficients a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its diffusion",
" coefficients.",
" The simplified scheme can be found in Kloeden, P.E. and Platen, E.:",
" Numerical Solution of Stochastic Differential Equations, pg. 458,",
" (Springer-Verlag 1992). The output consists of the variables YN[n],",
" DeltaWsM[n], and Delta[n]. YN[n] denotes the 1st order simplified",
" weak approximation to x[N] at the n-th step. DeltaWsM[n] denotes the",
" change in the M-dimensional noise process at the n-th step. Note here",
" that WsM[n] does not denote standard Wiener processes, but instead are",
" independent random variables (refer above text). Delta[n] denotes the",
" step size at the n-th step.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[wkeuler].",
" ",
"EXAMPLES: ",
" ",
">wkeuler([x[1],x[2]],[[r,0],[0,r]]);",
" ",
"table([",
" 1 = (Y1[n + 1] = Y1[n] + Y1[n] Delta[n] + r Delta Ws1[n])",
" 2 = (Y2[n + 1] = Y2[n] + Y2[n] Delta[n] + r Delta Ws2[n])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[LO], stochastic[LJ], stochastic[wktay2]"
):
stochastic[wktay2]:=
proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+a[i]*Delta[n]+1/2*L0(a[i],a,b)*Delta[n]^2+
sum('(op(j,op(i,b))+1/2*Delta[n]*(L0(op(j,op(i,b)),a,b)+
LJ(a[i],b,j)))*Delta*Ws.j[n]','j' = 1 .. nops(op(1,b)))+1/2*
sum('sum('LJ(op(j2,op(i,b)),b,j1)*(Delta^2*Ws.j1[n]*Ws.j2[n]+
V[j1,j2])','j1' = 1 .. nops(op(1,b)))',
'j2' = 1 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/wktay2` := TEXT(
"FUNCTION: stochastic[wktay2] - constructs stochastic Taylor schemes of",
"weak order 2.",
" ",
"CALLING SEQUENCE: wktay2([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call wktay2([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" simplified weak order 2 stochastic Taylor scheme for an N-dimensional",
" Stochastic Differential Equation (SDE) with M-dimensional noise which has",
" as its drift coefficients a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its",
" diffusion coefficients. The simplified scheme can be found in",
" Kloeden, P.E. and Platen, E.: Numerical Solution of Stochastic",
" Differential Equations, pg. 466, (Springer-Verlag 1992).",
" The output consists of the variables YN[n],",
" DeltaWsM[n], V[(j1,j2)], and Delta[n]. YN[n] denotes the 2nd order",
" simplified weak approximation to x[N] at the n-th step. DeltaWsM[n]",
" denotes the change in the M-dimensional noise process at the n-th step.",
" Note here that WsM[n] does not denote standard Wiener processes, but",
" instead are independent random variables (refer above text). V[(j1,j2)]",
" denotes independent two-point random variables (refer above text).",
" Delta[n] denotes the step size at the n-th step.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[wktay2].",
" ",
"EXAMPLES: ",
" ",
">wktay2([x[1],x[2]],[[r,0],[0,r]]);",
" ",
"table([",
" 2",
" 1 = (Y1[n + 1] = Y1[n] + Y1[n] Delta[n] + 1/2 Y1[n] Delta[n]",
" + (r + 1/2 Delta[n] r) Delta Ws1[n])",
" 2",
" 2 = (Y2[n + 1] = Y2[n] + Y2[n] Delta[n] + 1/2 Y2[n] Delta[n]",
" + (r + 1/2 Delta[n] r) Delta Ws2[n])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ], stochastic[wkeuler]"
):
stochastic[MLJ]:=
proc(X::algebraic,a::list(algebraic),b::list(list(algebraic)),j::integer)
local flag;
flag := 0;
if j = 0 then flag := L0(X,a,b) fi;
if flag = 0 then flag := sum('op(j,op(k,b))*diff(X,x[k])',
'k' = 1 .. nops(b)) fi;
RETURN(flag)
end:
`help/text/MLJ` := TEXT(
"FUNCTION: stochastic[MLJ] - applies one of the partial differential",
"operators, either L0 or LJ, to X. These operators can be found in Kloeden,",
"P.E., Platen, E.: Numerical Solution of Stochastic Differential Equations,",
"(Springer-Verlag, 1992).",
" ",
"CALLING SEQUENCE: MLJ(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]],j);",
" ",
"PARAMETERS: X - algebraic, given in the variables x[N] and t.",
" a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" j - integer.",
" ",
"SYNOPSIS: ",
"- The call MLJ(X,[a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]],j); computes the",
" application of either LJ or L0 to X. LJ and L0 are partial differential",
" operators, as described in Kloeden, P.E., Platen, E.: Numerical",
" Solution of Stochastic Differential Equations, pg 339,",
" (Springer-Verlag, 1992).",
" a1,..,aN are the drift coefficients of the N-dimensional Stochastic",
" Differential Equation (SDE).",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M, where M is the dimension of the Wiener process. j=1,..,M refers to",
" the ''current'' dimension of the Wiener process.",
" The output variables are consistent with the variables used as input.",
" ",
"- This routine is used by the routine stochastic[wktay3].",
" In general, it is not intended for use on its own.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[MLJ].",
" ",
"EXAMPLES: ",
" ",
">MLJ(x[2],[x[2],x[1]],[[r,s],[s,r]],2);",
" ",
" r ",
" ",
">MLJ(x[2],[x[2],x[1]],[[r,s],[s,r]],0);",
" ",
" x[1] ",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[LJ],",
" stochastic[SL0], stochastic[wktay3]."
):
stochastic[wktay3]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,soln;
for i to nops(a) do
soln[i] := Y.i[n+1] = Y.i[n]+a[i]*Delta[n]+
sum('op(j,op(i,b))*Delta*W.j[n]','j' = 1 .. nops(op(1,b)))+
sum('MLJ(a[i],a,b,j0)*I[j0,0]','j0' = 0 .. nops(op(1,b)))+
sum('sum('MLJ(op(j2,op(i,b)),a,b,j1)*I[j1,j2]',
'j2' = 1 .. nops(op(1,b)))','j1' = 0 .. nops(op(1,b)))+
sum('sum('MLJ(MLJ(a[i],a,b,k2),a,b,k1)*I[k1,k2,0]',
'k1' = 0 .. nops(op(1,b)))','k2' = 0 .. nops(op(1,b)))+sum(
'sum('sum('MLJ(MLJ(op(m3,op(i,b)),a,b,m2),a,b,m1)*I[m1,m2,m3]',
'm3' = 1 .. nops(op(1,b)))','m2' = 0 .. nops(op(1,b)))',
'm1' = 0 .. nops(op(1,b)));
for u to nops(a) do soln[i] := subs(x[u] = Y.u[n],soln[i]) od
od;
RETURN(eval(soln))
end:
`help/text/wktay3` := TEXT(
"FUNCTION: stochastic[wktay3] - constructs stochastic Taylor schemes of",
"weak order 3.",
" ",
"CALLING SEQUENCE: wktay3([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call wktay3([a1,..,aN],[[b11,..,b1M],..,[bN1,..,bNM]]); returns the",
" weak order 3 stochastic Taylor scheme for an N-dimensional",
" Stochastic Differential Equation (SDE) with M-dimensional noise which has",
" as its drift coefficients a1,..,aN and [b11,..,b1M],..,[bN1,..,bNM] as its",
" diffusion coefficients. The weak order 3 scheme can be found in",
" Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic Differential",
" Equations, pg. 468, (Springer-Verlag 1992). The output consists of the",
" variables YN[n], DeltaWM[n], I[(j1,j2)], I[(j1,j2,j3)], and Delta[n].",
" YN[n] denotes the 2nd order weak approximation to x[N] at the n-th step.",
" DeltaWM[n] denotes the change in the M-dimensional Wiener process",
" at the n-th step. I[(j1,j2)] and I[(j1,j2,j3)] denote multiple Ito",
" integrals (refer above text). Delta[n] denotes the step size at the n-th",
" step.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[wktay3].",
" ",
"EXAMPLES: ",
" ",
" ",
">wktay3([x[2],x[1]],[[r,s],[s,r]]);",
" ",
"table([",
" 1 = (Y1[n + 1] = Y1[n] + Y2[n] Delta[n] + r Delta W1[n] + s Delta W2[n]",
" + Y1[n] I[0, 0] + s I[1, 0] + r I[2, 0] + Y2[n] I[0, 0, 0]",
" + r I[1, 0, 0] + s I[2, 0, 0])",
" 2 = (Y2[n + 1] = Y2[n] + Y1[n] Delta[n] + s Delta W1[n] + r Delta W2[n]",
" + Y2[n] I[0, 0] + r I[1, 0] + s I[2, 0] + Y1[n] I[0, 0, 0]",
" + s I[1, 0, 0] + r I[2, 0, 0])",
" ])",
" ",
"SEE ALSO: stochastic, stochastic[L0], stochastic[MLJ],",
" stochastic[wkeuler], stochastic[wktay2]."
):
stochastic[colour]:=
proc(a::list(algebraic),b::list(algebraic))
local temp1,i;
for i to nops(a) do temp1[i] := dx[i][t] = a[i]*dt+b[i]*z[t]*dt od;
temp1[i] := dz[t] = -gamma*z[t]*dt+beta*dW[t];
RETURN(eval(temp1))
end:
`help/text/colour` := TEXT(
"FUNCTION: stochastic[colour] - converts Stochastic Differential Equations ",
"(SDEs) with white noise into their coloured noise form.",
" ",
"CALLING SEQUENCE: colour([a1,..,aN],[b1,..,bN]);",
" ",
"PARAMETERS: a1,..,aN - algebraics, given in the variables x[N] and t.",
" b1,..,bN - algebraics, given in the variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call colour([a1,..,aN],[b1,..,bN]); converts N-dimensional SDEs",
" with scalar white noise into their coloured noise form. Details of this",
" procedure, together with reasoning and numerical examples can be found",
" in Milshtein, G.N., Tret'yakov, M.V.: Numerical Solution of Differential",
" Equations with Coloured Noise, J. Stat. Physics, (to be released).",
" a1,..,aN denotes the drift coefficients of the N-dimensional SDE.",
" b1,..,bN denotes the diffusion coefficients of the SDE (scalar noise).",
" The output consists of the variables z, x[N], W, gamma, beta, and t.",
" z denotes the Ornstein-Uhlenbeck process, or exponentially correlated",
" coloured noise. x[N] denotes the state variable of the (N+1)-dimensional",
" SDE. W denotes a standard Wiener process. gamma and beta denote",
" parameters, obtainable from experimental data. t denotes time.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[colour].",
" ",
"EXAMPLES: ",
" ",
">colour([a(x[1],t)],[b(x[1],t)]);",
" ",
"table([",
" 1 = (dx[1][t] = a(x[1], t) dt + b(x[1], t) z[t] dt)",
" 2 = (dz[t] = - gamma z[t] dt + beta dW[t])",
" ])",
" ",
">colour([x[2],x[1]*(alpha-x[1]^2)-x[2]],[0,sigma*x[1]]);",
" ",
"table([",
" 1 = (dx[1][t] = x[2] dt)",
" 2",
" 2 = (dx[2][t] = (x[1] (alpha - x[1] ) - x[2]) dt + sigma x[1] z[t] dt)",
" 3 = (dz[t] = - gamma z[t] dt + beta dW[t])",
" ])",
" ",
"SEE ALSO: stochastic"
):
stochastic[comm1]:=proc()
local LJ1,LJ2,k,j1,j2,flag,p;
for p to nargs do
if type(args[p],list) <> true then
ERROR(`Expecting input to be an expression sequence of lists`) fi
od;
for k to nargs do
for j1 to nops(args[1]) do
for j2 to nops(args[1]) do
LJ1 := sum('op(j1,args[l])*diff(op(j2,args[k]),x[l])',
'l' = 1 .. nargs);
LJ2 := sum('op(j2,args[l])*diff(op(j1,args[k]),x[l])',
'l' = 1 .. nargs);
if LJ1 <> LJ2 then flag := 1 fi
od
od
od;
if flag = 1 then
RETURN("Commutative noise of the first kind doesn't exist for this system")
else RETURN("This system exhibits commutative noise of the first kind")
fi;
end:
`help/text/comm1` := TEXT(
"FUNCTION: stochastic[comm1] - informs the user if the diffusion matrix of a",
"Stochastic Differential Equation (SDE) has commutative noise of the first",
"kind.",
" ",
"CALLING SEQUENCE: comm1([b11,..,b1M],..,[bN1,..,bNM]);",
" ",
"PARAMETERS: [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call comm1([b11,..,b1M],..,[bN1,..,bNM]); returns the user with a",
" statement indicating whether or not the diffusion matrix of the SDE has",
" commutative noise of the first kind.",
" The commutativity condition of the first kind can be found in Kloeden,",
" P.E., Platen, E.: Numerical Solution of Stochastic Differential",
" Equations, pg 348, (Springer-Verlag, 1992).",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M, where N is the dimension of the SDE and M is the dimension of the",
" Wiener process.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[comm1]",
" ",
"EXAMPLES: ",
" ",
">comm1([1,0],[0,x[1]]);",
" ",
" Commutative noise of the first kind doesn't exist for this system",
" ",
"SEE ALSO: stochastic, stochastic[comm2]"
):
stochastic[comm2]:=
proc()
local LJ1LJ2,LJ2LJ1,k,p,j1,j2,j3,flag;
for p to nargs do
if type(args[p],list) <> true then
ERROR(`Expecting input to be an expression sequence of lists`) fi;
od;
for k to nargs do
for j1 to nops(args[1]) do
for j2 to nops(args[1]) do
for j3 to nops(args[1]) do
LJ1LJ2 := sum('op(j1,args[m])*diff(sum('op(j2,args[l])*
diff(op(j3,args[k]),x[l])','l' = 1 .. nargs),x[m])',
'm' = 1 .. nargs);
LJ2LJ1 := sum('op(j2,args[m])*diff(sum('op(j1,args[l])*
diff(op(j3,args[k]),x[l])','l' = 1 .. nargs),x[m])',
'm' = 1 .. nargs);
if LJ1LJ2 <> LJ2LJ1 then flag := 1 fi;
od;
od;
od;
od;
if flag = 1 then
RETURN("Commutative noise of the second kind doesn't exist for this system")
else RETURN("This system exhibits commutative noise of the second kind")
fi;
end:
`help/text/comm2` := TEXT(
"FUNCTION: stochastic[comm2] - informs the user if the diffusion matrix of a",
"Stochastic Differential Equation (SDE) has commutative noise of the second",
"kind.",
" ",
"CALLING SEQUENCE: comm2([b11,..,b1M],..,[bN1,..,bNM]);",
" ",
"PARAMETERS: [b11,..,b1M],..,[bN1,..,bNM] - lists of algebraics, given in",
" the variables x[N] and t.",
" ",
"SYNOPSIS: ",
"- The call comm2([b11,..,b1M],..,[bN1,..,bNM]); returns the user with a",
" statement indicating whether or not the diffusion matrix of the SDE has",
" commutative noise of the second kind.",
" The commutativity condition of the second kind can be found in",
" Kloeden, P.E., Platen, E.: Numerical Solution of Stochastic Differential",
" Equations, pg 355, (Springer-Verlag, 1992).",
" [b11,..,b1M],..,[bN1,..,bNM] denotes the diffusion matrix of dimension",
" N*M, where N is the dimension of the SDE and M is the dimension of the",
" Wiener process.",
" ",
"- This routine is part of the stochastic package and is usually loaded via",
" the call with(stochastic). It can also be invoked via the call",
" stochastic[comm2].",
" ",
"EXAMPLES: ",
" ",
">comm2([1,(x[2])^2*(x[1])^4],[0,(x[1])^2]);",
" ",
" Commutative noise of the second kind doesn't exist for this system",
" ",
"SEE ALSO: stochastic, stochastic[comm1]"
):
stochastic[itoformula]:=proc(U::list(algebraic),a::list(algebraic),
b::list(list(algebraic)))
local i,k,l0,lj,soln;
for i from 1 to nops(U) do
l0:=L0(U[i],a,b)*dt;
lj:=0;
for k from 1 to nops(op(1,b)) do
lj:=lj+LJ(U[i],b,k)*dW.k;
od;
soln[i]:=dX.i=l0 +lj;
od;
RETURN(eval(soln));
end:
stochastic[LFP]:=proc(P::algebraic,a::list(algebraic),b::list(list(algebraic)))
local part1,part2,part3;
part1 := diff(P,t);
part2 := sum('diff(a[k]*P,y[k])','k' = 1 .. nops(a));
part3 := 1/2*sum(
'sum('sum('diff(op(j,op(k,b))*op(j,op(l,b))*P,y[k],y[l])',
'j' = 1 .. nops(op(1,b)))','k' = 1 .. nops(a))','l' = 1 .. nops(a)
);
part1+part2-part3;
end:
stochastic[conv]:=proc(a::list(algebraic),b::list(list(algebraic)),c::algebraic)
local temp,i;
if c=ito then
for i from 1 to nops(a) do temp[i]:=op(i,a)-1/2*sum('sum('op(k,op(j,b))
*diff(op(k,op(i,b)),x[j])','k'=1..nops(op(1,b)))','j'=1..nops(a));
od;
elif c=strat then
for i from 1 to nops(a) do
temp[i]:=op(i,a)+1/2*sum('sum('op(k,op(j,b))
*diff(op(k,op(i,b)),x[j])','k'=1..nops(op(1,b)))','j'=1..nops(a));
od;
else
ERROR(`Must enter either ito or strat for the 3rd argument`)
fi;
RETURN(map(simplify,eval(temp)))
end:
stochastic[chainrule]:=proc(U::list(algebraic),a::list(algebraic),
b::list(list(algebraic)))
local i,k,l0,lj,soln;
for i from 1 to nops(U) do
l0:=SL0(U[i],a,b)*dt;
lj:=0;
for k from 1 to nops(b) do
lj:=lj+LJ(U[i],b,k)*odW.i;
od;
soln[i]:=dX.i=l0 +lj;
od;
RETURN(eval(soln));
end:
stochastic[linearize]:=proc(a::list(algebraic),b::list(list(algebraic)),c::list(algebraic))
local i,tempA,tempB,j,k,l;
tempA:=array(1..nops(a),1..nops(a));
for i from 1 to nops(a) do for j from 1 to nops(a)
do tempA[i,j]:=diff(op(i,a),x[j]);
od; od;
for i from 1 to nops(a) do for j from 1 to nops(a) do for l from 1 to nops(c)
do tempA[i,j]:=subs(x[l]=op(l,c),tempA[i,j]);
od; od; od;
for k from 1 to nops(op(1,b)) do tempB[k]:=array(1..nops(a),1..nops(a));
for i from 1 to nops(a) do for j from 1 to nops(a)
do tempB[k][i,j]:=diff(op(k,op(i,b)),x[j]);
od; od;
for i from 1 to nops(a) do for j from 1 to nops(a) do for l from 1 to nops(c)
do tempB[k][i,j]:=subs(x[l]=op(l,c),tempB[k][i,j]);
od; od; od;
od;
RETURN(A=map(simplify,convert(eval(tempA),matrix)),B=eval(tempB))
end:
stochastic[sphere]:=proc(a,b)
global q,q0,qk,h,hk;
local i,j,k,tempa,tempb,stempbs,N,tmp;
if type(a,array) then tmp:=convert(a,listlist);else tmp:=a; fi;
N:=nops(tmp);
hk:=evaln(hk); h:=evaln(h); qk:=evaln(qk);
q:=evaln(q); q0:=evaln(q0);tempa:=evaln(tempa);tempb:=evaln(tempb);
q0:=sum('sum('s[i]*a[i,j]','i'=1..N)*s[j]','j'=1..N);
for k from 1 to nops(b) do
qk[k]:=sum('sum('s[i]*b[k][i,j]','i'=1..N)*s[j]','j'=1..N);
od;
for k from 1 to nops(b) do
stempbs[k]:=sum('sum('s[i]*(b[k][i,j]+b[k][j,i])','i'=1..N)*s[j]','j'=1..N);
od;
q:=q0+ sum('0.5*stempbs[k]-qk[k]^2','k'=1..nops(b));
for i from 1 to N do
for j from 1 to N do
if (i=j) then tempa[i,i]:=a[i,i]-q0;
else tempa[i,j]:=a[i,j];
fi;
od;
od;
for i from 1 to N do
h[i]:=sum('tempa[i,j]','j'=1..N);
od;
for k from 1 to nops(b) do
for i from 1 to N do
for j from 1 to N do
if (i=j) then tempb[k][i,i]:=b[k][i,i]-qk[k];
else tempb[k][i,j]:=b[k][i,j];
fi;
od;
od;
od;
for k from 1 to nops(b) do
for i from 1 to N do
hk[k][i]:=sum('tempb[k][i,j]','j'=1..N);
od;
od;
end:
stochastic[position]:=proc(N,i,j)
global stelle;
stelle:=sum('N-k+1','k'=1..i-1)+j-i+1;
end:
stochastic[ap]:=proc(A)
local i,j,k,Atmp,N,counter;
global B1;
if type(A,array) then Atmp:=convert(A,listlist);else Atmp:=A; fi;
N:=nops(Atmp);
B1:=array(1..N*(N+1)/2,1..N*(N+1)/2);
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
B1[i,j]:=0;
od;
od;
counter:=0:
for i from 1 to N do
for j from i to N do
counter:=counter+1;
for k from 1 to N do
if (j<=k) then B1[counter,position(N,j,k)]:=A[i,k];
else B1[counter,position(N,k,j)]:=A[i,k];
fi;
od;
od;
od;
RETURN(evalm(B1));
end:
stochastic[pa]:=proc(A)
local i,j,k,Atmp,N,counter;
global B2;
if type(A,array) then Atmp:=convert(A,listlist);else Atmp:=A; fi;
N:=nops(Atmp);
B2:=array(1..N*(N+1)/2,1..N*(N+1)/2);
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
B2[i,j]:=0;
od;
od;
counter:=0:
for i from 1 to N do
for j from i to N do
counter:=counter+1;
for k from 1 to N do
if (i<=k) then B2[counter,position(N,i,k)]:=A[j,k];
else B2[counter,position(N,k,i)]:=A[j,k];
fi;
od;
od;
od;
RETURN(evalm(B2));
end:
stochastic[bpb]:=proc(B)
local i,j,k,l,Btmp,N,counter;
global B3;
if type(B,array) then Btmp:=convert(B,listlist);else Btmp:=B; fi;
N:=nops(Btmp);
B3:=array(1..N*(N+1)/2,1..N*(N+1)/2);
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
B3[i,j]:=0;
od;
od;
counter:=0:
for i from 1 to N do
for j from i to N do
counter:=counter+1;
for l from 1 to N do
for k from 1 to N do
if (k<=l) then
B3[counter,position(N,k,l)]:=B3[counter,position(N,k,l)]
+B[i,k]*B[j,l];
else
B3[counter,position(N,l,k)]:=B3[counter,position(N,l,k)]
+B[i,k]*B[j,l];
fi;
od;
od;
od;
od;
RETURN(evalm(B3));
end:
stochastic[momenteqn]:=proc(A,B)
local i,j,k,N,Btmp,Ctmp;
global New_A;
if type(A,array) then Btmp:=convert(A,listlist);else Btmp:=A; fi;
N:=nops(Btmp);
New_A:=array(1..N*(N+1)/2,1..N*(N+1)/2);
Ctmp:=array(1..N*(N+1)/2,1..N*(N+1)/2);
stochastic[ap](A);
stochastic[pa](A);
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
Ctmp[i,j]:=0;
od;
od;
for k from 1 to nops(B) do
stochastic[bpb](B[k]);
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
Ctmp[i,j]:=Ctmp[i,j]+B3[i,j];
od;
od;
od;
for i from 1 to N*(N+1)/2 do
for j from 1 to N*(N+1)/2 do
New_A[i,j]:=B1[i,j]+B2[i,j]+Ctmp[i,j];
od;
od;
RETURN(evalm(New_A));
end:
stochastic[pmatrix2pvector]:=proc(p)
local i,j,k,ptmp;
global pvector;
if type(p,array) then ptmp:=convert(p,listlist);else ptmp:=p; fi;
pvector:=array(1..nops(ptmp)*(nops(ptmp)+1)/2);
k:=0;
for i from 1 to nops(ptmp) do
if (i>1) then k:=k+(nops(ptmp)-i+2); fi;
for j from i to nops(ptmp) do
pvector[k+j-i+1]:=ptmp[i,j];
od;
od;
RETURN(eval(pvector));
end:
stochastic[pvector2pmatrix]:=proc(pvector)
local i,j,k,ptmp,N;
global p;
if type(pvector,array) then ptmp:=convert(pvector,list);else
ptmp:=pvector; fi;
N:=-1/2+sqrt(1/4+2*nops(ptmp));
p:=array(1..N,1..N);
k:=0;
for i from 1 to N do
if (i>1) then k:=k+(N-i+2); fi;
for j from i to N do
p[i,j]:=ptmp[k+j-i+1];
if (i<>j) then p[j,i]:=p[i,j]; fi;
od;
od;
RETURN(eval(p));
end:
stochastic[milcomm]:=proc(a::list(algebraic),b::list(list(algebraic)))
local u,i,l,soln;
for i to nops(a) do
soln[i]:=Y.i[n+1]=Y.i[n]+L0(x[i],a,b)*Delta[n]
+sum('LJ(x[i],b,j)*Delta.W.j[n]','j'=1..nops(op(1,b)))
+1/2*sum('sum('LJ(op(j2,op(i,b)),b,j1)*
(Delta.W.j1[n])*(Delta.W.j2[n])',
'j1'=1..nops(op(1,b)))','j2'=1..nops(op(1,b)));
for l to nops(op(1,b)) do
soln[i]:=subs((Delta.W.l[n])^2=((Delta.W.l[n])^2-Delta[n]),
soln[i]) od;
for u to nops(a) do
soln[i]:=subs(x[u]=Y.u[n],soln[i]);
od; od;
RETURN(eval(soln));
end: