From 501d31445590e4d27ce91d66c1070655ddadbb97 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 5 Apr 2018 23:38:05 -0700 Subject: [PATCH] Add non-diagonal SDE problems --- REQUIRE | 2 +- src/DiffEqProblemLibrary.jl | 8 +++-- src/sde_premade_problems.jl | 62 +++++++++++++++++++++++++++++++++---- 3 files changed, 62 insertions(+), 10 deletions(-) diff --git a/REQUIRE b/REQUIRE index 8427390..9c401c5 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,5 @@ julia 0.6 ParameterizedFunctions 2.0.0 DiffEqBase 3.0.3 -DiffEqPDEBase DiffEqOperators +DiffEqBiological diff --git a/src/DiffEqProblemLibrary.jl b/src/DiffEqProblemLibrary.jl index ade1f40..9d44c6c 100644 --- a/src/DiffEqProblemLibrary.jl +++ b/src/DiffEqProblemLibrary.jl @@ -2,7 +2,7 @@ __precompile__() module DiffEqProblemLibrary -using DiffEqBase, ParameterizedFunctions, DiffEqPDEBase, DiffEqOperators +using DiffEqBase, ParameterizedFunctions, DiffEqOperators, DiffEqBiological include("ode/ode_premade_problems.jl") include("dae_premade_problems.jl") @@ -19,8 +19,10 @@ export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear, SolverDiffEq, filament_prob #SDE Example Problems - export prob_sde_wave, prob_sde_linear, prob_sde_cubic, prob_sde_2Dlinear, prob_sde_lorenz, - prob_sde_2Dlinear, prob_sde_additive, prob_sde_additivesystem, oval2ModelExample + export prob_sde_wave, prob_sde_linear, prob_sde_cubic, prob_sde_2Dlinear, + prob_sde_lorenz, prob_sde_2Dlinear, prob_sde_additive, + prob_sde_additivesystem, oval2ModelExample, prob_sde_bistable, + prob_sde_bruss, prob_sde_oscilreact #SDE Stratonovich Example Problems export prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich diff --git a/src/sde_premade_problems.jl b/src/sde_premade_problems.jl index 108f866..281605f 100644 --- a/src/sde_premade_problems.jl +++ b/src/sde_premade_problems.jl @@ -91,11 +91,10 @@ u(u0,p,t,W_t)=\\arctan(\\frac{W_t}{10} + \\tan(u0)) """ prob_sde_wave = SDEProblem(f,σ,1.,(0.0,1.0)) -const sde_wave_α = 0.1 -const sde_wave_β = 0.05 -f = (u,p,t) -> sde_wave_β./sqrt.(1+t) - u./(2*(1+t)) -σ = (u,p,t) -> sde_wave_α*sde_wave_β./sqrt.(1+t) -(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0./sqrt.(1+t) + sde_wave_β*(t+sde_wave_α*W)./sqrt.(1+t) +f = (u,p,t) -> p[2]./sqrt.(1+t) - u./(2*(1+t)) +σ = (u,p,t) -> p[1]*p[2]./sqrt.(1+t) +p = (0.1,0.05) +(ff::typeof(f))(::Type{Val{:analytic}},u0,p,t,W) = u0./sqrt.(1+t) + p[2]*(t+p[1]*W)./sqrt.(1+t) """ Additive noise problem @@ -110,7 +109,7 @@ and initial condition u0=1.0 with α=0.1 and β=0.05, with solution u(u0,p,t,W_t)=\\frac{u0}{\\sqrt{1+t}} + \\frac{β(t+αW_t)}{\\sqrt{1+t}} ``` """ -prob_sde_additive = SDEProblem(f,σ,1.,(0.0,1.0)) +prob_sde_additive = SDEProblem(f,σ,1.,(0.0,1.0),p) const sde_wave_αvec = [0.1;0.1;0.1;0.1] const sde_wave_βvec = [0.5;0.25;0.125;0.1115] @@ -386,3 +385,54 @@ function generate_stiff_stoch_heat(D=1,k=1;N = 100) end SDEProblem(f,g,ones(N),(0.0,3.0),noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3))) end + +bistable_f(du,u,p,t) = du[1] = p[1] + p[2]*u[1]^4/(u[1]^4+11.9^4) - p[3]*u[1] +function bistable_g(du,u,p,t) + du[1,1] = .1*sqrt(p[1] + p[2]*u[1]^4/(u[1]^4+11.9^4)) + du[1,2] = -.1*sqrt(p[3]*u[1]) +end +p = (5.0,18.0,1.0) +""" +Bistable chemical reaction network with a semi-stable lower state. +""" +prob_sde_bistable = SDEProblem(bistable_f,bistable_g,[3.],(0.,300.),p, + noise_rate_prototype=zeros(1,2)) + +function bruss_f(du,u,p,t) + du[1] = p[1] + p[2]*u[1]*u[1]*u[2] - p[3]*u[1]-p[4]*u[1] + du[2] = p[3]*u[1] - p[2]*u[1]*u[1]*u[2] +end +function bruss_g(du,u,p,t) + du[1,1] = 0.15*sqrt(p[1]) + du[1,2] = 0.15*sqrt(p[2]*u[1]*u[1]*u[2]) + du[1,3] = -0.15*sqrt(p[3]*u[1]) + du[1,4] = -0.15*sqrt(p[4]*u[1]) + du[2,1] = 0 + du[2,2] = -0.15*sqrt(p[2]*u[1]*u[1]*u[2]) + du[2,3] = 0.15*sqrt(p[3]*2.5*u[1]) + du[2,4] = 0 +end +p = (1.0,1.0,2.5,1.0) +""" +Stochastic Brusselator +""" +prob_sde_bruss = SDEProblem(bruss_f,bruss_g,[3.,2.],(0.,100.),p, + noise_rate_prototype=zeros(2,4)) + +network = @reaction_network rnType begin + p1, (X,Y,Z) --> 0 + hill(X,p2,100.,-4), 0 --> Y + hill(Y,p3,100.,-4), 0 --> Z + hill(Z,p4,100.,-4), 0 --> X + hill(X,p5,100.,6), 0 --> R + hill(Y,p6,100.,4)*0.002, R --> 0 + p7, 0 --> S + R*p8, S --> SP + p9, SP + SP --> SP2 + p10, SP2 --> 0 +end p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 +p = (0.01,3.0,3.0,4.5,2.,15.,20.0,0.005,0.01,0.05) +""" +An oscillatory chemical reaction system +""" +prob_sde_oscilreact = SDEProblem(network, [200.,60.,120.,100.,50.,50.,50.], (0.,4000.),p)