Skip to content

Commit

Permalink
sqdt
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChrisRackauckas
ChrisRackauckas committed Jan 24, 2017
1 parent 862c67d commit 193f3c5
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 26 deletions.
20 changes: 10 additions & 10 deletions src/integrators/integrator_utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ end
local T::tType
local ΔW::randType
local ΔZ::randType
@unpack u,t,dt,T,timeseries,Ws,ts,adaptivealg,δ,discard_length,unstable_check,rands,sqdt,W,Z = integrator
@unpack u,t,dt,T,timeseries,Ws,ts,adaptivealg,δ,discard_length,unstable_check,rands,W,Z = integrator

integrator.opts.progress && (prog = Juno.ProgressBar(name=integrator.opts.progress_name))
if uType <: AbstractArray
Expand All @@ -42,8 +42,8 @@ end
iter = 0
max_stack_size = 0
max_stack_size2 = 0
ΔW = sqdt*next(rands) # Take one first
ΔZ = sqdt*next(rands) # Take one first
ΔW = integrator.sqdt*next(rands) # Take one first
ΔZ = integrator.sqdt*next(rands) # Take one first
end

@def sde_sritableaupreamble begin
Expand Down Expand Up @@ -141,19 +141,19 @@ end
if adaptivealg==:RSwM1
if !isempty(S₁)
dt,ΔW,ΔZ = pop!(S₁)
sqdt = sqrt(dt)
integrator.sqdt = sqrt(dt)
else # Stack is empty
c = min(integrator.opts.dtmax,dtnew)
dt = max(min(c,abs(T-t)),integrator.opts.dtmin)#abs to fix complex sqrt issue at end
#dt = min(c,abs(T-t))
sqdt = sqrt(dt)
ΔW = sqdt*next(rands)
ΔZ = sqdt*next(rands)
integrator.sqdt = sqrt(dt)
ΔW = integrator.sqdt*next(rands)
ΔZ = integrator.sqdt*next(rands)
end
elseif adaptivealg==:RSwM2 || adaptivealg==:RSwM3
c = min(integrator.opts.dtmax,dtnew)
dt = max(min(c,abs(T-t)),integrator.opts.dtmin) #abs to fix complex sqrt issue at end
sqdt = sqrt(dt)
integrator.sqdt = sqrt(dt)
if !(uType <: AbstractArray)
dttmp = 0.0; ΔW = 0.0; ΔZ = 0.0
else
Expand Down Expand Up @@ -258,7 +258,7 @@ end
else
W = W + ΔW
end
ΔW = sqdt*next(rands)
ΔW = integrator.sqdt*next(rands)
if !(typeof(integrator.alg) <: EM) || !(typeof(integrator.alg) <: RKMil)
if uType <: AbstractArray
for i in eachindex(u)
Expand All @@ -267,7 +267,7 @@ end
else
Z = Z + ΔZ
end
ΔZ = sqdt*next(rands)
ΔZ = integrator.sqdt*next(rands)
end
@sde_savevalues
end
Expand Down
8 changes: 4 additions & 4 deletions src/integrators/low_order.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,11 @@ function sde_solve{uType<:AbstractArray,uEltype,Nm1,N,tType,tTypeNoUnits,uEltype
integrator.g(t,u,L)
for i in eachindex(u)
K[i] = u[i] + dt*du1[i]
utilde[i] = K[i] + L[i]*sqdt
utilde[i] = K[i] + L[i]*integrator.sqdt
end
integrator.g(t,utilde,du2)
for i in eachindex(u)
u[i] = K[i]+L[i]*ΔW[i]+(du2[i]-L[i])./(2sqdt).*(ΔW[i].^2 - dt)
u[i] = K[i]+L[i]*ΔW[i]+(du2[i]-L[i])./(2integrator.sqdt).*(ΔW[i].^2 - dt)
end
@sde_loopfooter
end
Expand All @@ -52,8 +52,8 @@ function sde_solve{uType<:Number,uEltype,Nm1,N,tType,tTypeNoUnits,uEltypeNoUnits

K = u + dt.*integrator.f(t,u)
L = integrator.g(t,u)
utilde = K + L*sqdt
u = K+L*ΔW+(integrator.g(t,utilde)-integrator.g(t,u))/(2sqdt)*(ΔW^2 - dt)
utilde = K + L*integrator.sqdt
u = K+L*ΔW+(integrator.g(t,utilde)-integrator.g(t,u))/(2integrator.sqdt)*(ΔW^2 - dt)

@sde_loopfooter
end
Expand Down
24 changes: 12 additions & 12 deletions src/integrators/sri.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ function sde_solve{algType<:SRI,uType<:AbstractArray,uEltype,Nm1,N,tType,tTypeNo
@sde_loopheader

for i in eachindex(u)
chi1[i] = .5*(ΔW[i].^2 - dt)/sqdt #I_(1,1)/sqrt(h)
chi1[i] = .5*(ΔW[i].^2 - dt)/integrator.sqdt #I_(1,1)/sqrt(h)
chi2[i] = .5*(ΔW[i] + ΔZ[i]/sqrt(3)) #I_(1,0)/h
chi3[i] = 1/6 * (ΔW[i].^3 - 3*ΔW[i]*dt)/dt #I_(1,1,1)/h
end
Expand All @@ -47,7 +47,7 @@ function sde_solve{algType<:SRI,uType<:AbstractArray,uEltype,Nm1,N,tType,tTypeNo
end
end
H0[i] = u + A0temp*dt + B0temp.*chi2
H1[i] = u + A1temp*dt + B1temp*sqdt
H1[i] = u + A1temp*dt + B1temp*integrator.sqdt
end
atemp[:]=zero(uEltype)
btemp[:]=zero(uEltype)
Expand Down Expand Up @@ -116,7 +116,7 @@ function sde_solve{uType<:AbstractArray,uEltype,Nm1,N,tType,tTypeNoUnits,uEltype
@sde_loopheader

for i in eachindex(u)
chi1[i] = (ΔW[i].^2 - dt)/2sqdt #I_(1,1)/sqrt(h)
chi1[i] = (ΔW[i].^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
chi2[i] = (ΔW[i] + ΔZ[i]/sqrt(3))/2 #I_(1,0)/h
chi3[i] = (ΔW[i].^3 - 3ΔW[i]*dt)/6dt #I_(1,1,1)/h
end
Expand All @@ -130,13 +130,13 @@ function sde_solve{uType<:AbstractArray,uEltype,Nm1,N,tType,tTypeNoUnits,uEltype
fH01o4[i] = fH01[i]/4
g₁o2[i] = g₁[i]/2
H0[i] = u[i] + 3*(fH01o4[i] + chi2[i]*g₁o2[i])
H11[i] = u[i] + fH01o4[i] + sqdt*g₁o2[i]
H12[i] = u[i] + fH01[i] - sqdt*g₁[i]
H11[i] = u[i] + fH01o4[i] + integrator.sqdt*g₁o2[i]
H12[i] = u[i] + fH01[i] - integrator.sqdt*g₁[i]
end
integrator.g(t+dto4,H11,g₂)
integrator.g(t+dt,H12,g₃)
for i in eachindex(u)
H13[i] = u[i] + fH01o4[i] + sqdt*(-5g₁[i] + 3g₂[i] + g₃[i]/2)
H13[i] = u[i] + fH01o4[i] + integrator.sqdt*(-5g₁[i] + 3g₂[i] + g₃[i]/2)
end

integrator.g(t+dto4,H13,g₄)
Expand Down Expand Up @@ -185,7 +185,7 @@ function sde_solve{uType<:Number,uEltype,Nm1,N,tType,tTypeNoUnits,uEltypeNoUnits
@inbounds while t<T
@sde_loopheader

chi1 = (ΔW.^2 - dt)/2sqdt #I_(1,1)/sqrt(h)
chi1 = (ΔW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
chi2 = (ΔW + ΔZ/sqrt(3))/2 #I_(1,0)/h
chi3 = (ΔW.^3 - 3ΔW*dt)/6dt #I_(1,1,1)/h
fH01 = dt*integrator.f(t,u)
Expand All @@ -195,11 +195,11 @@ function sde_solve{uType<:Number,uEltype,Nm1,N,tType,tTypeNoUnits,uEltypeNoUnits
dto4 = dt/4
g₁o2 = g₁/2
H0 = u + 3*(fH01o4 + chi2.*g₁o2)
H11 = u + fH01o4 + sqdt*g₁o2
H12 = u + fH01 - sqdt*g₁
H11 = u + fH01o4 + integrator.sqdt*g₁o2
H12 = u + fH01 - integrator.sqdt*g₁
g₂ = integrator.g(t+dto4,H11)
g₃ = integrator.g(t+dt,H12)
H13 = u + fH01o4 + sqdt*(-5g₁ + 3g₂ + g₃/2)
H13 = u + fH01o4 + integrator.sqdt*(-5g₁ + 3g₂ + g₃/2)


g₄ = integrator.g(t+dto4,H13)
Expand Down Expand Up @@ -241,7 +241,7 @@ function sde_solve{algType<:SRI,uType<:Number,uEltype,Nm1,N,tType,tTypeNoUnits,u
@inbounds while t<T
@sde_loopheader

chi1 = .5*(ΔW.^2 - dt)/sqdt #I_(1,1)/sqrt(h)
chi1 = .5*(ΔW.^2 - dt)/integrator.sqdt #I_(1,1)/sqrt(h)
chi2 = .5*(ΔW + ΔZ/sqrt(3)) #I_(1,0)/h
chi3 = 1/6 * (ΔW.^3 - 3*ΔW*dt)/dt #I_(1,1,1)/h

Expand All @@ -259,7 +259,7 @@ function sde_solve{algType<:SRI,uType<:Number,uEltype,Nm1,N,tType,tTypeNoUnits,u
B1temp += B₁[i,j]*integrator.g(t + c₁[j]*dt,H1[j])
end
H0[i] = u + A0temp*dt + B0temp.*chi2
H1[i] = u + A1temp*dt + B1temp*sqdt
H1[i] = u + A1temp*dt + B1temp*integrator.sqdt
end
atemp = zero(u)
btemp = zero(u)
Expand Down

0 comments on commit 193f3c5

Please sign in to comment.