From e4abed817a10d0e22c2cb5e3fbe9305caf69c144 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Wed, 26 Jun 2024 11:40:17 -0700 Subject: [PATCH 1/8] W2Ito1 scheme --- src/StochasticDiffEq.jl | 2 +- src/alg_utils.jl | 3 + src/algorithms.jl | 15 +- src/caches/srk_weak_caches.jl | 164 ++++++++++++++++++ src/perform_step/srk_weak.jl | 267 +++++++++++++++++++++++++++++ src/solve.jl | 12 ++ test/ode_convergence_regression.jl | 3 + test/runtests.jl | 1 + test/weak_convergence/W2Ito1.jl | 228 ++++++++++++++++++++++++ 9 files changed, 693 insertions(+), 2 deletions(-) create mode 100644 test/weak_convergence/W2Ito1.jl diff --git a/src/StochasticDiffEq.jl b/src/StochasticDiffEq.jl index 61cff9659..3e98fb5ca 100644 --- a/src/StochasticDiffEq.jl +++ b/src/StochasticDiffEq.jl @@ -151,7 +151,7 @@ end SROCK1, SROCK2, SROCKEM, SKSROCK, TangXiaoSROCK2, KomBurSROCK2, SROCKC2, WangLi3SMil_A, WangLi3SMil_B, WangLi3SMil_C, WangLi3SMil_D, WangLi3SMil_E, WangLi3SMil_F, AutoSOSRI2, AutoSOSRA2, - DRI1, DRI1NM, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM, + DRI1, DRI1NM, RI1, RI3, RI5, RI6, RDI1WM, RDI2WM, RDI3WM, RDI4WM, W2Ito1, RS1, RS2, PL1WM, PL1WMA, NON, COM, NON2 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 0d49bd483..f97824e6b 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -101,6 +101,7 @@ alg_order(alg::RDI1WM) = 1 // 1 alg_order(alg::RDI2WM) = 1 // 1 alg_order(alg::RDI3WM) = 1 // 1 alg_order(alg::RDI4WM) = 1 // 1 +alg_order(alg::W2Ito1) = 1 // 1 alg_order(alg::RS1) = 1 // 1 alg_order(alg::RS2) = 1 // 1 @@ -181,6 +182,7 @@ alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RDI1WM) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RDI2WM) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RDI3WM) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RDI4WM) = true +alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::W2Ito1) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RS1) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::RS2) = true alg_compatible(prob::DiffEqBase.AbstractSDEProblem, alg::PL1WM) = true @@ -254,6 +256,7 @@ alg_needs_extra_process(alg::RDI1WM) = true alg_needs_extra_process(alg::RDI2WM) = true alg_needs_extra_process(alg::RDI3WM) = true alg_needs_extra_process(alg::RDI4WM) = true +alg_needs_extra_process(alg::W2Ito1) = true alg_needs_extra_process(alg::RS1) = true alg_needs_extra_process(alg::RS2) = true alg_needs_extra_process(alg::PL1WM) = true diff --git a/src/algorithms.jl b/src/algorithms.jl index b38800595..743ace4d8 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -93,7 +93,7 @@ Defaults to solving the Ito problem, but RKMilCommute(interpretation=:Stratonovi Uses a 1.5/2.0 error estimate for adaptive time stepping. Default: ii_approx=IICommutative() does not approximate the Levy area. """ -struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm +struct RKMilCommute{T} <: StochasticDiffEqAdaptiveAlgorithm interpretation::Symbol ii_approx::T end @@ -483,6 +483,19 @@ Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise. """ struct RDI4WM <: StochasticDiffEqAdaptiveAlgorithm end + +""" +Tang, X., & Xiao, A., Efficient weak second-order stochastic Runge–Kutta methods +for Itô stochastic differential equations, +BIT Numerical Mathematics, 57, 241-260 (2017) +DOI: 10.1007/s10543-016-0618-9 + +W2Ito1: High Weak Order Method +Adaptive step weak order 2.0 for Ito SDEs (deterministic order 3). +Can handle diagonal, non-diagonal, non-commuting, and scalar additive noise. +""" +struct W2Ito1 <: StochasticDiffEqAdaptiveAlgorithm end + # Stratonovich sense """ diff --git a/src/caches/srk_weak_caches.jl b/src/caches/srk_weak_caches.jl index 9334e0459..23328f843 100644 --- a/src/caches/srk_weak_caches.jl +++ b/src/caches/srk_weak_caches.jl @@ -2318,3 +2318,167 @@ function alg_cache(alg::SMEB,prob,u,ΔW,ΔZ,p,rate_prototype, SIESMECache(u,uprev,W2,W3,tab,k0,k1,g0,g1,g2,tmpu) end + + + + + +# Tang & Xiao: DOI 10.1007/s10543-016-0618-9 W2Ito1 and W2Ito2 methods + +struct W2Ito1ConstantCache{T} <: StochasticDiffEqConstantCache + # hard-coded version + a021::T + a031::T + a032::T + + a121::T + a131::T + #a132::T + + #a221::T + #a231::T + #a232::T + + b021::T + b031::T + #b032::T + + b121::T + b131::T + #b132::T + + b221::T + #b222::T + #b223::T + #b231::T + #b232::T + #b233::T + + α1::T + α2::T + α3::T + + beta01::T + beta02::T + beta03::T + + beta11::T + #beta12::T + beta13::T + + #quantile(Normal(),1/6) + NORMAL_ONESIX_QUANTILE::T +end + + + +function W2Ito1ConstantCache(::Type{T}, ::Type{T2}) where {T,T2} + + a021 = convert(T, 1 // 2) + a031 = convert(T, -1) + a032 = convert(T, 2) + + a121 = convert(T, 1 // 4) + a131 = convert(T, 1 // 4) + + b021 = convert(T, (6 - sqrt(6)) / 10) + b031 = convert(T, (3 + 2 * sqrt(6)) / 5) + + b121 = convert(T, 1 // 2) + b131 = convert(T, -1 // 2) + + b221 = convert(T, 1) + + α1 = convert(T, 1 // 6) + α2 = convert(T, 2 // 3) + α3 = convert(T, 1 // 6) + + beta01 = convert(T, -1) + beta02 = convert(T, 1) + beta03 = convert(T, 1) + + beta11 = convert(T, 2) + beta13 = convert(T, -2) + + NORMAL_ONESIX_QUANTILE = convert(T, -0.9674215661017014) + + W2Ito1ConstantCache(a021, a031, a032, a121, a131, b021, b031, b121, b131, b221, α1, α2, α3, beta01, beta02, beta03, beta11, beta13, NORMAL_ONESIX_QUANTILE) +end + + +function alg_cache(alg::W2Ito1, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} + W2Ito1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) +end + +@cache struct W2Ito1Cache{uType,randType,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache + u::uType + uprev::uType + uhat::uType + + _dW::randType + _dZ::randType + chi1::randType + + tab::tabType + + g1::rateNoiseType + g2::rateNoiseType + g3::rateNoiseType + + k1::rateType + k2::rateType + k3::rateType + + H02::uType + H03::uType + H12::Vector{uType} + H13::Vector{uType} + + tmp1::possibleRateType + tmpg::rateNoiseType + + tmp::uType + resids::uType + +end + +function alg_cache(alg::W2Ito1, prob, u, ΔW, ΔZ, p, rate_prototype, + noise_rate_prototype, jump_rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} + if ΔW isa Union{SArray,Number} + _dW = copy(ΔW) + _dZ = zeros(eltype(ΔW), 2) + chi1 = copy(ΔW) + else + _dW = zero(ΔW) + _dZ = zeros(eltype(ΔW), 2) + chi1 = zero(ΔW) + end + m = length(ΔW) + tab = W2Ito1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits)) + g1 = zero(noise_rate_prototype) + g2 = zero(noise_rate_prototype) + g3 = zero(noise_rate_prototype) + k1 = zero(rate_prototype) + k2 = zero(rate_prototype) + k3 = zero(rate_prototype) + + H02 = zero(u) + H03 = zero(u) + H12 = Vector{typeof(u)}() + H13 = Vector{typeof(u)}() + + for k = 1:m + push!(H12, zero(u)) + push!(H13, zero(u)) + end + + tmp1 = zero(rate_prototype) + tmpg = zero(noise_rate_prototype) + + uhat = copy(uprev) + tmp = zero(u) + resids = zero(u) + + W2Ito1Cache(u, uprev, uhat, _dW, _dZ, chi1, tab, g1, g2, g3, k1, k2, k3, H02, H03, H12, H13, tmp1, tmpg, tmp, resids) +end diff --git a/src/perform_step/srk_weak.jl b/src/perform_step/srk_weak.jl index 125739817..2bf78ea1f 100644 --- a/src/perform_step/srk_weak.jl +++ b/src/perform_step/srk_weak.jl @@ -2343,3 +2343,270 @@ end @. u = u + γ1*g0*W.dW + (λ1*W.dW + λ2*integrator.sqdt + λ3*W2)*g1 + (µ1*W.dW + µ2*integrator.sqdt + µ3*W2)*g2 end end + + +# W2Ito1 due to Tang and Xiao + +function W2Ito_I2(_dW, xi, eta2, k, l) + if k < l + return 1 // 2 * (_dW[l] - eta2 * _dW[l]) + elseif k > l + return 1 // 2 * (_dW[l] + eta2 * _dW[l]) + else + return 1 // 2 * (_dW[k]^2 / xi - xi) + end +end + +@muladd function perform_step!(integrator, cache::W2Ito1ConstantCache) + @unpack a021, a031, a032, a121, a131, b021, b031, b121, b131, b221, α1, α2, α3, beta01, beta02, beta03, beta11, beta13, NORMAL_ONESIX_QUANTILE = cache + @unpack t, dt, uprev, u, W, p, f = integrator + + # define three-point distributed random variables + dW_scaled = W.dW / sqrt(dt) + sq3dt = sqrt(3 * dt) + _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), dW_scaled) + # define two-point distributed random variables + _dZ = map(x -> calc_twopoint_random(one(integrator.dt), x), W.dZ) + xi = integrator.sqdt * _dZ[1] + + m = length(W.dW) + + chi1 = map(x -> (x^2 / xi - xi) / 2, _dW) # diagonal of Ihat2 + + # compute stage values + k1 = integrator.f(uprev, p, t) + g1 = integrator.g(uprev, p, t) + + # H_1^(0) (stage 1) + # H01 = uprev + # H_2^(0) (stage 2) + if !is_diagonal_noise(integrator.sol.prob) || W.dW isa Number + H02 = uprev + a021 * k1 * dt + b021 * g1 * _dW + else + H02 = uprev + a021 * k1 * dt + b021 * g1 .* _dW + end + + # H_3^(0) (stage 3) + k2 = integrator.f(H02, p, t) + H03 = uprev + a032 * k2 * dt + a031 * k1 * dt + + if !is_diagonal_noise(integrator.sol.prob) || W.dW isa Number + H03 += b031 * g1 * _dW + else + H03 += b031 * g1 .* _dW + end + + k3 = integrator.f(H03, p, t) + + + # H_1^(k), (stage 1) + # H11 = uprev + # H_i^(k) (stage 2 and 3) + if W.dW isa Number + H12 = uprev + a121 * k1 * dt + b121 * g1 * xi + H13 = uprev + a131 * k1 * dt + b131 * g1 * xi + else + H12 = [uprev .+ a121 * k1 * dt for k = 1:m] + H13 = [uprev .+ a131 * k1 * dt for k = 1:m] + for k = 1:m + if is_diagonal_noise(integrator.sol.prob) + tmp = zero(integrator.u) + tmp[k] = g1[k] + H12[k] += b121 * tmp * xi + H13[k] += b131 * tmp * xi + for l = 1:m + if l!=k + tmp = zero(integrator.u) + tmp[l] = g1[l] + H12[k] += b221 * tmp * W2Ito_I2(_dW, xi, _dZ[2], k, l) + end + end + + else + H12[k] += b121 * g1[:, k] * xi + H13[k] += b131 * g1[:, k] * xi + + for l = 1:m + if l!=k + H12[k] += b221 * g1[:, l] * W2Ito_I2(_dW, xi, _dZ[2], k, l) + end + end + end + end + end + + if W.dW isa Number + g2 = integrator.g(H12, p, t) + g3 = integrator.g(H13, p, t) + else + if is_diagonal_noise(integrator.sol.prob) + g2 = [integrator.g(H12[k], p, t)[k] for k = 1:m] + g3 = [integrator.g(H13[k], p, t)[k] for k = 1:m] + else + g2 = hcat([integrator.g(H12[k], p, t)[:, k] for k = 1:m] ...) + g3 = hcat([integrator.g(H13[k], p, t)[:, k] for k = 1:m] ...) + end + end + + # add stages together Eq. (3) + u = uprev + α1 * k1 * dt + α2 * k2 * dt + α3 * k3 * dt + + # add noise + if W.dW isa Number + u += g1 * (_dW * beta01 + chi1 * beta11) + g2 * (_dW * beta02) + g3 * (_dW * beta03 + chi1 * beta13) + else + if is_diagonal_noise(integrator.sol.prob) + u += g1 .* (_dW * beta01 + chi1 * beta11) + g2 .* (_dW * beta02) + g3 .* (_dW * beta03 + chi1 * beta13) + else + # non-diag noise + u += g1 * (_dW * beta01 + chi1 * beta11) + g2 * (_dW * beta02) + g3 * (_dW * beta03 + chi1 * beta13) + end + end + + if integrator.opts.adaptive + + # schemes with lower convergence order + # check against EM + uhat = uprev + k1 * dt + + if is_diagonal_noise(integrator.sol.prob) + uhat += g1 .* _dW + else + uhat += g1 * _dW + end + + resids = calculate_residuals(u - uhat, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(resids, t) + end + + integrator.u = u + +end + + +@muladd function perform_step!(integrator, cache::W2Ito1Cache) + @unpack t, dt, uprev, u, W, p, f = integrator + @unpack _dW, _dZ, chi1, tab, g1, g2, g3, k1, k2, k3, H02, H03, H12, H13, tmp1, tmpg, uhat, tmp, resids = cache + @unpack a021, a031, a032, a121, a131, b021, b031, b121, b131, b221, α1, α2, α3, beta01, beta02, beta03, beta11, beta13, NORMAL_ONESIX_QUANTILE = cache.tab + + m = length(W.dW) + sq3dt = sqrt(3 * dt) + + if W.dW isa Union{SArray,Number} + # define three-point distributed random variables + _dW = map(x -> calc_threepoint_random(sq3dt, NORMAL_ONESIX_QUANTILE, x), W.dW / sqrt(dt)) + # define two-point distributed random variables + _dZ = map(x -> calc_twopoint_random(one(integrator.dt), x), W.dZ) + xi = integrator.sqdt * _dZ[1] + chi1 = map(x -> (x^2 / xi - xi) / 2, _dW) # diagonal of Ihat2 + else + # define three-point distributed random variables + sqrtdt = sqrt(dt) + @.. chi1 = W.dW / sqrtdt + calc_threepoint_random!(_dW, sq3dt, NORMAL_ONESIX_QUANTILE, chi1) + calc_twopoint_random!(_dZ, one(integrator.sqdt), W.dZ) + xi = integrator.sqdt * _dZ[1] + map!(x -> (x^2 / xi - xi) / 2, chi1, _dW) + end + + # compute stage values + integrator.f(k1, uprev, p, t) + integrator.g(g1, uprev, p, t) + + # H_i^(0), stage 1 + # H01 = uprev + # H_i^(0), stage 2 + if is_diagonal_noise(integrator.sol.prob) + @.. H02 = uprev + a021 * k1 * dt + b021 * g1 * _dW + else + mul!(tmp1, g1, _dW) + @.. H02 = uprev + dt * a021 * k1 + b021 * tmp1 + end + + integrator.f(k2, H02, p, t) + if is_diagonal_noise(integrator.sol.prob) + @.. H03 = uprev + a032 * k2 * dt + a031 * k1 * dt + b031 * g1 * _dW + else + @.. H03 = uprev + a032 * k2 * dt + a031 * k1 * dt + b031 * tmp1 + end + + integrator.f(k3, H03, p, t) + + # H_i^(k), stages + # H11 = uprev + for k = 1:m + if is_diagonal_noise(integrator.sol.prob) + fill!(tmpg, zero(eltype(integrator.u))) + tmpg[k] = g1[k] + @.. H12[k] = uprev + a121 * k1 * dt + b121 * tmpg * xi + @.. H13[k] = uprev + a131 * k1 * dt + b131 * tmpg * xi + for l = 1:m + if l!=k + fill!(tmpg, zero(eltype(integrator.u))) + tmpg[l] = g1[l] + WikJ = W2Ito_I2(_dW, xi, _dZ[2], k, l) + @.. H12[k] = H12[k] + b221 * tmpg * WikJ + end + end + integrator.g(tmpg, H12[k], p, t) + g2[k] = tmpg[k] + integrator.g(tmpg, H13[k], p, t) + g3[k] = tmpg[k] + else + g1k = @view g1[:, k] + @.. H12[k] = uprev + a121 * k1 * dt + b121 * g1k * xi + @.. H13[k] = uprev + a131 * k1 * dt + b131 * g1k * xi + for l = 1:m + if l!=k + WikJ = W2Ito_I2(_dW, xi, _dZ[2], k, l) + g1l = @view g1[:, l] + @.. H12[k] = H12[k] + b221 * g1l * WikJ + end + end + integrator.g(tmpg, H12[k], p, t) + tmpgk = @view tmpg[:, k] + g2k = @view g2[:, k] + copyto!(g2k, tmpgk) + integrator.g(tmpg, H13[k], p, t) + tmpgk = @view tmpg[:, k] + g3k = @view g3[:, k] + copyto!(g3k, tmpgk) + end + end + + # add stages together Eq. (3) + @.. u = uprev + α1 * k1 * dt + α2 * k2 * dt + α3 * k3 * dt + + # add noise + if W.dW isa Number || is_diagonal_noise(integrator.sol.prob) + @.. u = u + g1 * (_dW * beta01 + chi1 * beta11) + g2 * (_dW * beta02) + g3 * (_dW * beta03 + chi1 * beta13) + else + # non-diag noise + mul!(tmp1, g1, (_dW * beta01 + chi1 * beta11)) + @.. u = u + tmp1 + mul!(tmp1, g2, (_dW * beta02)) + @.. u = u + tmp1 + mul!(tmp1, g3, (_dW * beta03 + chi1 * beta13)) + @.. u = u + tmp1 + end + + + if integrator.opts.adaptive + + # check against EM + @.. uhat = uprev + k1 * dt + + if is_diagonal_noise(integrator.sol.prob) + @.. uhat = uhat + g1 .* _dW + else + mul!(tmp1, g1, _dW) + @.. uhat = uhat + tmp1 + end + @.. tmp = u - uhat + + calculate_residuals!(resids, tmp, uprev, uhat, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t) + + integrator.EEst = integrator.opts.internalnorm(resids, t) + end +end diff --git a/src/solve.jl b/src/solve.jl index d06354c81..c940e479a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -345,6 +345,12 @@ function DiffEqBase.__init( W = WienerProcess!(t,rand_prototype,rand_prototype2, save_everystep=save_noise, rng = Xorshifts.Xoroshiro128Plus(_seed)) + elseif alg isa W2Ito1 + m = 2 + rand_prototype2 = zeros(eltype(rand_prototype), Int(m)) + W = WienerProcess!(t, rand_prototype, rand_prototype2, + save_everystep=save_noise, + rng=Xorshifts.Xoroshiro128Plus(_seed)) else W = WienerProcess!(t,rand_prototype,rand_prototype, save_everystep=save_noise, @@ -393,6 +399,12 @@ function DiffEqBase.__init( W = WienerProcess(t,rand_prototype,rand_prototype2, save_everystep=save_noise, rng = Xorshifts.Xoroshiro128Plus(_seed)) + elseif alg isa W2Ito1 + m = 2 + rand_prototype2 = zeros(eltype(rand_prototype), Int(m)) + W = WienerProcess(t, rand_prototype, rand_prototype2, + save_everystep=save_noise, + rng=Xorshifts.Xoroshiro128Plus(_seed)) else W = WienerProcess(t,rand_prototype,rand_prototype, save_everystep=save_noise, diff --git a/test/ode_convergence_regression.jl b/test/ode_convergence_regression.jl index 01d65ea68..80c0f1dde 100644 --- a/test/ode_convergence_regression.jl +++ b/test/ode_convergence_regression.jl @@ -51,6 +51,9 @@ sim2 = test_convergence(dts,prob, RDI3WM(),trajectories=20) sim2 = test_convergence(dts,prob, RDI4WM(),trajectories=20) @test abs(sim2.𝒪est[:l∞]-3) <.1 +sim2 = test_convergence(dts, prob, W2Ito1(), trajectories=20) +@test abs(sim2.𝒪est[:l∞] - 3) < 0.1 + sim2 = test_convergence(dts,prob, RS1(),trajectories=20) @test abs(sim2.𝒪est[:l∞]-2) <.1 sim2 = test_convergence(dts,prob, RS2(),trajectories=20) diff --git a/test/runtests.jl b/test/runtests.jl index abe68baa4..fe5fb8911 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -86,6 +86,7 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") end if !is_APPVEYOR && GROUP == "WeakConvergence3" + @time @safetestset "Tang & Xiao weak SRK Tests" begin include("weak_convergence/W2Ito1.jl") end @time @safetestset "Roessler weak SRK (non-diagonal) Tests" begin include("weak_convergence/srk_weak_final_non_diagonal.jl") end end diff --git a/test/weak_convergence/W2Ito1.jl b/test/weak_convergence/W2Ito1.jl new file mode 100644 index 000000000..be332d60d --- /dev/null +++ b/test/weak_convergence/W2Ito1.jl @@ -0,0 +1,228 @@ +""" + Tests for W2Ito1 from Tang, X., & Xiao, A. (2017). Efficient weak second-order stochastic + Runge–Kutta methods for Itô stochastic differential equations. BIT Numerical Mathematics, + 57, 241-260. +""" + +import Statistics # for mean values of trajectories +import LinearAlgebra # for the normn +using StochasticDiffEq +using Test +using Random +using DiffEqDevTools +seed = 103473 +function prob_func(prob, i, repeat) + remake(prob, seed=seeds[i]) +end + +""" + Test OOP +""" + +@info "Scalar noise" + +numtraj = Int(1e6) # in the paper they use 1e9 +u₀ = 0.0 +f(u, p, t) = 1 // 2 * u + sqrt(u^2 + 1) +g(u, p, t) = sqrt(u^2 + 1) +dts = 1 .// 2 .^ (4:-1:1) +tspan = (0.0, 2.0) # 2.0 in paper + + +h1(z) = z^3 - 6 * z^2 + 8 * z +#analytical_sol(t) = E(f(X(t))) = E(h1(arsinh(X(t))) = t^3-3*t^2+2*t +#analytical_sol(2) = 0 and analytical_sol(1)=0 + +Random.seed!(seed) +seeds = rand(UInt, numtraj) + +prob = SDEProblem(f, g, u₀, tspan) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h1(asinh(sol[end])), false), + prob_func=prob_func +) + + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=0.0 +) +@test abs(sim.𝒪est[:weak_final] - 2) < 0.3 +println("W2Ito1:", sim.𝒪est[:weak_final]) + +@info "Diagonal noise" + +u₀ = [0.1, 0.1] +function f2(u, p, t) + [3 // 2 * u[1], 3 // 2 * u[2]] +end +function g2(u, p, t) + [1 // 10 * u[1], 1 // 10 * u[2]] +end +dts = 1 .// 2 .^ (3:-1:0) +tspan = (0.0, 1.0) + +h2(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) if h3(z) = z^2 ) + +prob = SDEProblem(f2, g2, u₀, tspan) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h2(sol[end][1]), false), + prob_func=prob_func +) + +numtraj = Int(1e5) +Random.seed!(seed) +seeds = rand(UInt, numtraj) + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=1 // 100 * exp(301 // 100) +) +@test -(sim.𝒪est[:weak_final] - 2) < 0.3 # order is 2.4 +println("W2Ito1:", sim.𝒪est[:weak_final]) + + +@info "Non-commutative noise" + +u₀ = [1.0, 1.0] +function f3(u, p, t) + return [-273 // 512 * u[1], -1 // 160 * u[1] - (-785 // 512 + sqrt(2) / 8) * u[2]] +end +function g3(u, p, t) + [1 // 4 * u[1] 1 // 16 * u[1] + (1 - 2 * sqrt(2)) / 4 * u[1] 1 // 10 * u[1] + 1 // 16 * u[2]] +end +dts = 1 .// 2 .^ (3:-1:0) +tspan = (0.0, 3.0) + +h3(z) = z^2 # but apply it only to u[1] + +prob = SDEProblem(f3, g3, u₀, tspan, noise_rate_prototype=zeros(2, 2)) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h3(sol[end][1]), false), + prob_func=prob_func +) + +numtraj = Int(1e5) +seed = 100 +Random.seed!(seed) +seeds = rand(UInt, numtraj) + + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=exp(-3.0) +) +@test abs(sim.𝒪est[:weak_final] - 2) < 0.3 + +println("W2Ito1:", sim.𝒪est[:weak_final]) + +# """ +# Test IIP +# """ + +@info "Scalar noise" + +u₀ = [0.0] +f1!(du, u, p, t) = @.(du = 1 // 2 * u + sqrt(u^2 + 1)) +g1!(du, u, p, t) = @.(du = sqrt(u^2 + 1)) +dts = 1 .// 2 .^ (4:-1:1) +tspan = (0.0, 2.0) + +h1(z) = z^3 - 6 * z^2 + 8 * z + +prob = SDEProblem(f1!, g1!, u₀, tspan) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h1(asinh(sol[end][1])), false), + prob_func=prob_func +) + +numtraj = Int(1e6) +Random.seed!(seed) +seeds = rand(UInt, numtraj) + + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=0.0 +) +@test abs(sim.𝒪est[:weak_final] - 2) < 0.3 +println("W2Ito1:", sim.𝒪est[:weak_final]) + +@info "Diagonal noise" + +u₀ = [0.1, 0.1] +function f2!(du, u, p, t) + du[1] = 3 // 2 * u[1] + du[2] = 3 // 2 * u[2] +end +function g2!(du, u, p, t) + du[1] = 1 // 10 * u[1] + du[2] = 1 // 10 * u[2] +end +dts = 1 .// 2 .^ (3:-1:0) +tspan = (0.0, 1.0) + +h2(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) if h3(z) = z^2 ) + +prob = SDEProblem(f2!, g2!, u₀, tspan) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h2(sol[end][1]), false), + prob_func=prob_func +) + +numtraj = Int(1e5) +Random.seed!(seed) +seeds = rand(UInt, numtraj) + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=1 // 100 * exp(301 // 100) +) +@test -(sim.𝒪est[:weak_final] - 2) < 0.3 # order is 2.4 +println("W2Ito1:", sim.𝒪est[:weak_final]) + + +@info "Non-commutative noise" + +u₀ = [1.0, 1.0] +function f3!(du, u, p, t) + du[1] = -273 // 512 * u[1] + du[2] = -1 // 160 * u[1] - (-785 // 512 + sqrt(2) / 8) * u[2] +end +function g3!(du, u, p, t) + du[1, 1] = 1//4*u[1] + du[1, 2] = 1//16*u[1] + du[2, 1] = (1-2*sqrt(2))/4*u[1] + du[2, 2] = 1//10*u[1]+1//16*u[2] +end +dts = 1 .// 2 .^ (3:-1:0) +tspan = (0.0, 3.0) + +h3(z) = z^2 # but apply it only to u[1] + +prob = SDEProblem(f3!, g3!, u₀, tspan, noise_rate_prototype=zeros(2, 2)) +ensemble_prob = EnsembleProblem(prob; + output_func=(sol, i) -> (h3(sol[end][1]), false), + prob_func=prob_func +) + +numtraj = Int(1e5) +seed = 100 +Random.seed!(seed) +seeds = rand(UInt, numtraj) + + +sim = test_convergence(dts, ensemble_prob, W2Ito1(), + save_everystep=false, trajectories=numtraj, save_start=false, adaptive=false, + weak_timeseries_errors=false, weak_dense_errors=false, + expected_value=exp(-3.0) +) +@test abs(sim.𝒪est[:weak_final] - 2) < 0.3 + +println("W2Ito1:", sim.𝒪est[:weak_final]) From 6bcf3406e9cc2d55854552c090d142df3817fad8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 30 Oct 2024 15:13:44 -0100 Subject: [PATCH 2/8] Update pipeline.yml --- .buildkite/pipeline.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 3fb870af3..7f6ddafc0 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -12,7 +12,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "1" + version: "lts" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -32,7 +32,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "1" + version: "lts" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -53,7 +53,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "1" + version: "lts" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -68,7 +68,7 @@ steps: - label: "WeakAdaptiveGPU" plugins: - JuliaCI/julia#v1: - version: "1" + version: "lts" - JuliaCI/julia-test#v1: coverage: false agents: From 3ac4d71e67de776f8a8fe883c16b65f7eecb3883 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 31 Oct 2024 17:24:03 -0100 Subject: [PATCH 3/8] Update pipeline.yml --- .buildkite/pipeline.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7f6ddafc0..7fdf77425 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -81,5 +81,6 @@ steps: if: build.message !~ /\[skip tests\]/ env: + debug_plugin: "true" JULIA_PKG_SERVER: "" # it often struggles with our large artifacts # SECRET_CODECOV_TOKEN: "..." From 4db5be2389b7df2581c03ab8b7a4e2583485cdee Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 31 Oct 2024 17:32:54 -0100 Subject: [PATCH 4/8] Update pipeline.yml --- .buildkite/pipeline.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 7fdf77425..3ec921558 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,6 +20,7 @@ steps: os: "linux" queue: "juliaecosystem" arch: "x86_64" + debug_plugin: "true" timeout_in_minutes: 240 # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ @@ -81,6 +82,5 @@ steps: if: build.message !~ /\[skip tests\]/ env: - debug_plugin: "true" JULIA_PKG_SERVER: "" # it often struggles with our large artifacts # SECRET_CODECOV_TOKEN: "..." From 7ff23ee1fdd743a1316fd6b866ab2d424f8fb89c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 31 Oct 2024 17:35:11 -0100 Subject: [PATCH 5/8] Update pipeline.yml --- .buildkite/pipeline.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 3ec921558..0730d9470 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -12,7 +12,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "lts" + version: "1.10" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -33,7 +33,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "lts" + version: "1.10" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -54,7 +54,7 @@ steps: GROUP: "{{matrix}}" plugins: - JuliaCI/julia#v1: - version: "lts" + version: "1.10" - JuliaCI/julia-test#v1: coverage: false julia_args: "--threads=auto" @@ -69,7 +69,7 @@ steps: - label: "WeakAdaptiveGPU" plugins: - JuliaCI/julia#v1: - version: "lts" + version: "1.10" - JuliaCI/julia-test#v1: coverage: false agents: From ed21ce868a5c68bab1b20f2cebde80eb67df5a4c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 3 Nov 2024 18:09:12 -0100 Subject: [PATCH 6/8] Update pipeline.yml --- .buildkite/pipeline.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 0730d9470..b11e35cb9 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -20,7 +20,6 @@ steps: os: "linux" queue: "juliaecosystem" arch: "x86_64" - debug_plugin: "true" timeout_in_minutes: 240 # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ From 0dd11d36e5b577d24545412f6726c80290ddef9e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sun, 3 Nov 2024 20:15:41 -0100 Subject: [PATCH 7/8] Update test/weak_convergence/W2Ito1.jl --- test/weak_convergence/W2Ito1.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/weak_convergence/W2Ito1.jl b/test/weak_convergence/W2Ito1.jl index be332d60d..3c98a843f 100644 --- a/test/weak_convergence/W2Ito1.jl +++ b/test/weak_convergence/W2Ito1.jl @@ -208,7 +208,7 @@ h3(z) = z^2 # but apply it only to u[1] prob = SDEProblem(f3!, g3!, u₀, tspan, noise_rate_prototype=zeros(2, 2)) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h3(sol[end][1]), false), + output_func=(sol, i) -> (h3(sol.u[end][1]), false), prob_func=prob_func ) From 1ef01f8b69c6d31e9b2644d5d5fbdcc39524bd9f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Mon, 4 Nov 2024 00:25:14 -0100 Subject: [PATCH 8/8] depwarn fixes --- test/adaptive/sde_weak_adaptive.jl | 2 +- test/complex_tests.jl | 4 ++-- test/gpu/sde_weak_adaptive_gpu.jl | 2 +- test/rode_linear_tests.jl | 20 +++++++++---------- test/split_tests.jl | 2 +- test/tdir_tests.jl | 2 +- test/weak_convergence/PL1WM.jl | 12 +++++------ test/weak_convergence/SIE_SME.jl | 6 +++--- test/weak_convergence/W2Ito1.jl | 10 +++++----- .../srk_weak_diagonal_final.jl | 2 +- test/weak_convergence/srk_weak_final.jl | 4 ++-- .../srk_weak_final_non_diagonal.jl | 2 +- test/weak_convergence/weak_strat.jl | 6 +++--- .../weak_strat_non_diagonal.jl | 2 +- 14 files changed, 38 insertions(+), 38 deletions(-) diff --git a/test/adaptive/sde_weak_adaptive.jl b/test/adaptive/sde_weak_adaptive.jl index 62a1d02d7..cfc88a30b 100644 --- a/test/adaptive/sde_weak_adaptive.jl +++ b/test/adaptive/sde_weak_adaptive.jl @@ -23,7 +23,7 @@ function reduction(u,batch,I) end function output_func(sol,i) - #h1(asinh(sol[end][1])),false + #h1(asinh(sol.u[end][1])),false h1.(asinh.(sol)),false end diff --git a/test/complex_tests.jl b/test/complex_tests.jl index 9e4143de4..1b14e2c7c 100644 --- a/test/complex_tests.jl +++ b/test/complex_tests.jl @@ -15,7 +15,7 @@ implicit_noautodiff = [SKenCarp(autodiff=false), ImplicitEulerHeun(autodiff=fals for alg in Iterators.flatten((explicit, implicit_noautodiff)) sol = solve(prob, alg) - @test eltype(sol[end]) == ComplexF64 + @test eltype(sol.u[end]) == ComplexF64 end # currently broken @@ -34,7 +34,7 @@ end for alg in Iterators.flatten((explicit, implicit_noautodiff)) sol = solve(prob, alg) - @test eltype(sol[end]) == ComplexF64 + @test eltype(sol.u[end]) == ComplexF64 end # currently broken diff --git a/test/gpu/sde_weak_adaptive_gpu.jl b/test/gpu/sde_weak_adaptive_gpu.jl index cfa415e55..bfe0e46d0 100644 --- a/test/gpu/sde_weak_adaptive_gpu.jl +++ b/test/gpu/sde_weak_adaptive_gpu.jl @@ -25,7 +25,7 @@ function reduction(u,batch,I) end function output_func(sol,i) - #h1(asinh(sol[end][1])),false + #h1(asinh(sol.u[end][1])),false h1.(asinh.(sol)),false end diff --git a/test/rode_linear_tests.jl b/test/rode_linear_tests.jl index 735ea134e..eb79a89a2 100644 --- a/test/rode_linear_tests.jl +++ b/test/rode_linear_tests.jl @@ -8,9 +8,9 @@ prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) -@test abs(sol[end] - sol2[end]) < 0.1 * abs(sol[end]) +@test abs(sol.u[end] - sol2[end]) < 0.1 * abs(sol.u[end]) sol3 = solve(prob,RandomTamedEM(),dt=1/100) -@test abs(sol[end] - sol3[end]) < 0.1 * abs(sol[end]) +@test abs(sol.u[end] - sol3[end]) < 0.1 * abs(sol.u[end]) f(du,u,p,t,W) = (du.=1.01u.+0.87u.*W) u0 = ones(4) @@ -18,9 +18,9 @@ prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) -@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) +@test sum(abs,sol.u[end]-sol2[end]) < 0.1 * sum(abs, sol.u[end]) sol3 = solve(prob,RandomEM(),dt=1/100) -@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol[end]) +@test sum(abs,sol.u[end]-sol3[end]) < 0.1 * sum(abs, sol.u[end]) f(u,p,t,W) = 2u*sin(W) u0 = 1.00 @@ -29,9 +29,9 @@ prob = RODEProblem{false}(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) prob = RODEProblem{false}(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) -@test abs(sol[end]-sol2[end]) < 0.1 * abs(sol[end]) +@test abs(sol.u[end]-sol2[end]) < 0.1 * abs(sol.u[end]) sol3 = solve(prob,RandomTamedEM(),dt=1/100) -@test abs(sol[end]-sol3[end]) < 0.1 * abs(sol[end]) +@test abs(sol.u[end]-sol3[end]) < 0.1 * abs(sol.u[end]) function f(du,u,p,t,W) du[1] = 0.2u[1]*sin(W[1] - W[2]) @@ -43,9 +43,9 @@ prob = RODEProblem(f,u0,tspan) sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) -@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) +@test sum(abs,sol.u[end]-sol2[end]) < 0.1 * sum(abs, sol.u[end]) sol3 = solve(prob,RandomTamedEM(),dt=1/100) -@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol[end]) +@test sum(abs,sol.u[end]-sol3[end]) < 0.1 * sum(abs, sol.u[end]) function f(du,u,p,t,W) du[1] = -0.2W[3]*u[1]*sin(W[1] - W[2]) @@ -57,6 +57,6 @@ prob = RODEProblem(f,u0,tspan,rand_prototype=zeros(3)) sol = solve(prob,RandomEM(),dt=1/100, save_noise=true) prob = RODEProblem(f,u0,tspan, noise=NoiseWrapper(sol.W)) sol2 = solve(prob,RandomHeun(),dt=1/100) -@test sum(abs,sol[end]-sol2[end]) < 0.1 * sum(abs, sol[end]) +@test sum(abs,sol.u[end]-sol2[end]) < 0.1 * sum(abs, sol.u[end]) sol3 = solve(prob,RandomTamedEM(),dt=1/100) -@test sum(abs,sol[end]-sol3[end]) < 0.1 * sum(abs, sol3[end]) +@test sum(abs,sol.u[end]-sol3[end]) < 0.1 * sum(abs, sol3[end]) diff --git a/test/split_tests.jl b/test/split_tests.jl index efa4704b7..36c3f1bad 100644 --- a/test/split_tests.jl +++ b/test/split_tests.jl @@ -25,7 +25,7 @@ prob = SDEProblem{false}(f,σ,u0,(0.0,1.0),noise = NoiseWrapper(sol.W)) sol2 = solve(prob,EM(),dt=1/10) -@test sol[end][:] ≈ sol2[end][:] +@test sol.u[end][:] ≈ sol2[end][:] ################################################################################ diff --git a/test/tdir_tests.jl b/test/tdir_tests.jl index c51c09227..8ccfe06c9 100644 --- a/test/tdir_tests.jl +++ b/test/tdir_tests.jl @@ -26,6 +26,6 @@ prob = SDEProblem(f,g,u₀,(0.0,1.0)) sol =solve(prob,EulerHeun(),dt=0.01,save_noise=true) _sol = deepcopy(sol) # to make sure the plot is correct W3 = NoiseGrid(reverse!(_sol.t),reverse!(_sol.W)) -prob3 = SDEProblem(f,g,sol[end],(1.0,0.0),noise=W3) +prob3 = SDEProblem(f,g,sol.u[end],(1.0,0.0),noise=W3) sol2 = solve(prob3,EulerHeun(),dt=0.01) @test sol.u ≈ reverse!(sol2.u) atol=1e-1 diff --git a/test/weak_convergence/PL1WM.jl b/test/weak_convergence/PL1WM.jl index d9e073707..dea45493f 100644 --- a/test/weak_convergence/PL1WM.jl +++ b/test/weak_convergence/PL1WM.jl @@ -38,7 +38,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(f,g,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end]),false), + output_func = (sol,i) -> (h1(sol.u[end]),false), prob_func = prob_func ) @@ -67,7 +67,7 @@ p = [3//2,1//100] prob = SDEProblem(f1!,g1!,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end][1]),false), + output_func = (sol,i) -> (h1(sol.u[end][1]),false), prob_func = prob_func ) @@ -124,7 +124,7 @@ h2(z) = z prob = SDEProblem(f2!,g2!,u₀,tspan,p,noise_rate_prototype=zeros(4,4)) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h2(sol[end][1]),false), + output_func = (sol,i) -> (h2(sol.u[end][1]),false), prob_func = prob_func ) @@ -164,7 +164,7 @@ h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f3!,g3!,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h3(sol[end][1]),false), + output_func = (sol,i) -> (h3(sol.u[end][1]),false), prob_func = prob_func ) @@ -200,7 +200,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(fadd,gadd,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (sol[end],false), + output_func = (sol,i) -> (sol.u[end],false), prob_func = prob_func ) @@ -233,7 +233,7 @@ gadd!(du,u,p,t) = @.(du = p[2]) prob = SDEProblem(fadd!,gadd!,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (sol[end][1],false), + output_func = (sol,i) -> (sol.u[end][1],false), prob_func = prob_func ) diff --git a/test/weak_convergence/SIE_SME.jl b/test/weak_convergence/SIE_SME.jl index 88afc2e78..0ea320575 100644 --- a/test/weak_convergence/SIE_SME.jl +++ b/test/weak_convergence/SIE_SME.jl @@ -38,7 +38,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(f,g,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end]),false), + output_func = (sol,i) -> (h1(sol.u[end]),false), prob_func = prob_func ) @@ -90,7 +90,7 @@ p = [3//2,1//100] prob = SDEProblem(f1!,g1!,u₀,tspan,p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end][1]),false), + output_func = (sol,i) -> (h1(sol.u[end][1]),false), prob_func = prob_func ) @@ -155,7 +155,7 @@ h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f3!,g3!,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h3(sol[end][1]),false), + output_func = (sol,i) -> (h3(sol.u[end][1]),false), prob_func = prob_func ) diff --git a/test/weak_convergence/W2Ito1.jl b/test/weak_convergence/W2Ito1.jl index 3c98a843f..232ea22e0 100644 --- a/test/weak_convergence/W2Ito1.jl +++ b/test/weak_convergence/W2Ito1.jl @@ -38,7 +38,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(f, g, u₀, tspan) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h1(asinh(sol[end])), false), + output_func=(sol, i) -> (h1(asinh(sol.u[end])), false), prob_func=prob_func ) @@ -67,7 +67,7 @@ h2(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f2, g2, u₀, tspan) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h2(sol[end][1]), false), + output_func=(sol, i) -> (h2(sol.u[end][1]), false), prob_func=prob_func ) @@ -101,7 +101,7 @@ h3(z) = z^2 # but apply it only to u[1] prob = SDEProblem(f3, g3, u₀, tspan, noise_rate_prototype=zeros(2, 2)) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h3(sol[end][1]), false), + output_func=(sol, i) -> (h3(sol.u[end][1]), false), prob_func=prob_func ) @@ -136,7 +136,7 @@ h1(z) = z^3 - 6 * z^2 + 8 * z prob = SDEProblem(f1!, g1!, u₀, tspan) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h1(asinh(sol[end][1])), false), + output_func=(sol, i) -> (h1(asinh(sol.u[end][1])), false), prob_func=prob_func ) @@ -171,7 +171,7 @@ h2(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f2!, g2!, u₀, tspan) ensemble_prob = EnsembleProblem(prob; - output_func=(sol, i) -> (h2(sol[end][1]), false), + output_func=(sol, i) -> (h2(sol.u[end][1]), false), prob_func=prob_func ) diff --git a/test/weak_convergence/srk_weak_diagonal_final.jl b/test/weak_convergence/srk_weak_diagonal_final.jl index 6a0558fc0..8d631b0ff 100644 --- a/test/weak_convergence/srk_weak_diagonal_final.jl +++ b/test/weak_convergence/srk_weak_diagonal_final.jl @@ -37,7 +37,7 @@ h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f3!,g3!,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h3(sol[end][1]),false), + output_func = (sol,i) -> (h3(sol.u[end][1]),false), prob_func = prob_func ) diff --git a/test/weak_convergence/srk_weak_final.jl b/test/weak_convergence/srk_weak_final.jl index 682067185..fde70631b 100644 --- a/test/weak_convergence/srk_weak_final.jl +++ b/test/weak_convergence/srk_weak_final.jl @@ -38,7 +38,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(f,g,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(asinh(sol[end])),false), + output_func = (sol,i) -> (h1(asinh(sol.u[end])),false), prob_func = prob_func ) @@ -157,7 +157,7 @@ h1(z) = z^3-6*z^2+8*z prob = SDEProblem(f1!,g1!,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(asinh(sol[end][1])),false), + output_func = (sol,i) -> (h1(asinh(sol.u[end][1])),false), prob_func = prob_func ) diff --git a/test/weak_convergence/srk_weak_final_non_diagonal.jl b/test/weak_convergence/srk_weak_final_non_diagonal.jl index cc5363ee4..ec36bd71e 100644 --- a/test/weak_convergence/srk_weak_final_non_diagonal.jl +++ b/test/weak_convergence/srk_weak_final_non_diagonal.jl @@ -42,7 +42,7 @@ h2(z) = z^2 # but apply it only to u[1] prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2)) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h2(sol[end][1]),false), + output_func = (sol,i) -> (h2(sol.u[end][1]),false), prob_func = prob_func ) diff --git a/test/weak_convergence/weak_strat.jl b/test/weak_convergence/weak_strat.jl index f5fc32d82..625d5f8de 100644 --- a/test/weak_convergence/weak_strat.jl +++ b/test/weak_convergence/weak_strat.jl @@ -41,7 +41,7 @@ seeds = rand(UInt, numtraj) prob = SDEProblem(f,g,u₀,tspan, p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end]),false), + output_func = (sol,i) -> (h1(sol.u[end]),false), prob_func = prob_func ) @@ -115,7 +115,7 @@ dts = 1 .//2 .^(5:-1:0) prob = SDEProblem(f!,g!,[u₀],tspan, p) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h1(sol[end][1]),false), + output_func = (sol,i) -> (h1(sol.u[end][1]),false), prob_func = prob_func ) @@ -203,7 +203,7 @@ h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and == 1//100**exp(301//100*t) prob = SDEProblem(f3!,g3!,u₀,tspan) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h3(sol[end][1]),false), + output_func = (sol,i) -> (h3(sol.u[end][1]),false), prob_func = prob_func ) diff --git a/test/weak_convergence/weak_strat_non_diagonal.jl b/test/weak_convergence/weak_strat_non_diagonal.jl index 5a4c42473..8f935f924 100644 --- a/test/weak_convergence/weak_strat_non_diagonal.jl +++ b/test/weak_convergence/weak_strat_non_diagonal.jl @@ -43,7 +43,7 @@ h2(z) = z^2 # E(x_i) = 1/10 exp(1/2t) or E(x_1* x_2) = 1/100 exp(2t) prob = SDEProblem(f2!,g2!,u₀,tspan,noise_rate_prototype=zeros(2,2)) ensemble_prob = EnsembleProblem(prob; - output_func = (sol,i) -> (h2(sol[end][1]),false), + output_func = (sol,i) -> (h2(sol.u[end][1]),false), prob_func = prob_func )