Merge pull request #62 from frankschae/noise_grid_problem

noise wrapper alternative

Author: ChrisRackauckas

src/convergence.jl

@@ -61,28 +61,46 @@ function analyticless_test_convergence(dts::AbstractArray,
                           alg,test_dt;trajectories=100,
                           save_everystep=true,timeseries_steps=1,
                           timeseries_errors=save_everystep,adaptive=false,
-                          weak_timeseries_errors=false,weak_dense_errors=false,kwargs...)
+                          weak_timeseries_errors=false,weak_dense_errors=false,use_noise_grid=true,kwargs...)
   _solutions = []
   tmp_solutions = Array{Any}(undef,trajectories,length(dts))
   for j in 1:trajectories
     @info "Monte Carlo iteration: $j/$trajectories"
     t = prob.tspan[1]:test_dt:prob.tspan[2]
-    if prob.noise_rate_prototype === nothing
-      brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
-      brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+    if use_noise_grid
+      if prob.noise_rate_prototype === nothing
+        brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+        brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
+      else
+        brownian_values = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
+        brownian_values2 = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
+      end
+      np = NoiseGrid(t,brownian_values,brownian_values2)
+      _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,
+                         noise=np,
+                         noise_rate_prototype=prob.noise_rate_prototype);
+      true_sol = solve(_prob,alg;adaptive=adaptive,dt=test_dt);
+
+      for i in 1:length(dts)
+        sol = solve(_prob,alg;dt=dts[i],adaptive=adaptive);
+        err_sol = appxtrue(sol,true_sol)
+        tmp_solutions[j,i] = err_sol
+      end
     else
-      brownian_values = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
-      brownian_values2 = cumsum([[zeros(size(prob.noise_rate_prototype,2))];[sqrt(test_dt)*randn(size(prob.noise_rate_prototype,2)) for i in 1:length(t)-1]])
-    end
-    np = NoiseGrid(t,brownian_values,brownian_values2)
-    _prob = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,
-                       noise=np,
-                       noise_rate_prototype=prob.noise_rate_prototype);
-    true_sol = solve(_prob,alg;adaptive=adaptive,dt=test_dt);
-    for i in 1:length(dts)
-      sol = solve(_prob,alg;dt=dts[i],adaptive=adaptive);
-      err_sol = appxtrue(sol,true_sol)
-      tmp_solutions[j,i] = err_sol
+      # using NoiseWrapper doesn't lead to constant true_sol
+      true_sol = solve(prob,alg;adaptive=adaptive,dt=test_dt, save_noise=true)
+      _sol = deepcopy(true_sol)
+      W1 = NoiseWrapper(_sol.W)
+
+      _prob = remake(prob,u0=prob.u0,p=prob.p,tspan=prob.tspan,noise=W1,noise_rate_prototype = prob.noise_rate_prototype)
+
+      for i in 1:length(dts)
+        W1 = NoiseWrapper(_sol.W)
+        _prob = remake(prob,u0=prob.u0,p=prob.p,tspan=prob.tspan,noise=W1,noise_rate_prototype = prob.noise_rate_prototype)
+        sol = solve(_prob,alg;dt=dts[i],adaptive=adaptive);
+        err_sol = appxtrue(sol,true_sol)
+        tmp_solutions[j,i] = err_sol
+      end
     end
   end
   _solutions = [EnsembleSolution(tmp_solutions[:,i],0.0,true) for i in 1:length(dts)]

test/analyticless_convergence_tests.jl

@@ -14,21 +14,28 @@ test_setup = Dict(:alg=>Vern9(),:reltol=>1e-25,:abstol=>1e-25)
 sim1 = analyticless_test_convergence(dts,prob,Tsit5(),test_setup)
 sim2 = analyticless_test_convergence(dts,prob,Vern9(),test_setup)
 
-@test sim1.𝒪est[:final]-5 < 0.2
-@test sim2.𝒪est[:final]-9 < 0.2
+@test abs(sim1.𝒪est[:final]-5) < 0.2
+@test abs(sim2.𝒪est[:final]-9) < 0.2
 
 function f2(du,u,p,t)
   du .= 1.01u
 end
 function g2(du,u,p,t)
   du .= 1.01u
 end
-prob = SDEProblem(f2,g2,[1.0;1.0],(0.0,10.0))
+prob = SDEProblem(f2,g2,[1.0;1.0],(0.0,1.0))
 
 using StochasticDiffEq
 
-dts = (1/2).^(6:-1:3)
+dts = (1/2).^(7:-1:3)
 test_dt = 1/2^8
 Random.seed!(100)
-sim1 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=400)
-@test sim1.𝒪est[:final]-1.5 < 0.3
+sim1 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100)
+@test abs(sim1.𝒪est[:final]-1.5) < 0.4
+@show sim1.𝒪est[:final]
+
+dts = (1/2).^(7:-1:4)
+test_dt = 1/2^8
+sim2 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=100, use_noise_grid=false)
+@test abs(sim2.𝒪est[:final]-1.5) < 0.3
+@show sim2.𝒪est[:final]