diff --git a/lib/BVProblemLibrary/src/BVProblemLibrary.jl b/lib/BVProblemLibrary/src/BVProblemLibrary.jl index b66f198..cb3c0c7 100644 --- a/lib/BVProblemLibrary/src/BVProblemLibrary.jl +++ b/lib/BVProblemLibrary/src/BVProblemLibrary.jl @@ -5,6 +5,365 @@ using DiffEqBase, Markdown include("linear.jl") include("nonlinear.jl") +################### flat_moon ############################ +function flat_moon_function!(du, u, p, t) + g = 1.62 + A = 3 * 1.62 + du[1] = u[3] + du[2] = u[4] + du[3] = A * cos(u[5]) + du[4] = (A * sin(u[5]) - g) + du[5] = -u[6] * cos(u[5]) + du[6] = (u[6])^2 * sin(u[5]) + du[7] = 0 +end + +function flat_moon_bc!(res, u, p, t) + Vc = 1627 + h = 185.2 + res[1] = u[1][1] + res[2] = u[1][2] + res[3] = u[1][3] + res[4] = u[1][4] + res[5] = u[end][2] - h + res[6] = u[end][3] - Vc + res[7] = u[end][4] +end +tspan = (0, 700) +@doc raw""" + flat_moon + +This test problem is about the optimal-time launching of a satellite into orbit from flat moon without atmospheric drag. + +Given by + +```math +\frac{dz_1}{dt}=z_3t_f +``` +```math +\frac{dz_2}{dt}=z_4t_f +``` +```math +\frac{dz_3}{dt}=A\cos(z_5)t_f +``` +```math +\frac{dz_4}{dt}=(A\sin(z_5)-g)t_f +``` +```math +\frac{dz_5}{dt}=-z_6\cos(z_5)t_F +``` +```math +\frac{dz_6}{dt}=z_6^2\sin(z_5)t_f +``` +```math +\frac{dz_7}{dt}=0 +``` + +with boundary condition + +```math +z_1(0)=0 +``` +```math +z_2(0)=0 +``` +```math +z_3(0)=0 +``` +```math +z_4(0)=0 +``` +```math +z_5(1)=h +``` +```math +z_6(1)=V_c +``` +```math +z_7(1)=0 +``` + +# Solution + +No analytical solution + +# References + +[Reference](https://archimede.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=534) +""" +flat_moon = BVProblem(flat_moon_function!, flat_moon_bc!, [0, 0, 0, 0, 0, 0, 0], tspan) + +################### flat_earth ############################ +function flat_earth_function!(du, u, p, t) + Vc = sqrt(398600.4 / (6378.14 + 300)) * 1000 + h = 300000 + g = 9.80665 + thrust2Weight = 3 + acc = Thrust2Weight * g + + du[1] = (u[3] * (Vc / h)) + du[2] = (u[4] * (Vc / h)) + du[3] = (acc * (1 / (abs(Vc) * sqrt(1 + u[6]^2)))) + du[4] = (acc * (u[6] / (abs(Vc) * sqrt(1 + u[6]^2))) - (g / Vc)) + du[5] = 0 + du[6] = (-u[5] * (Vc / h)) + du[7] = 0 +end + +function flat_earth_bc!(res, u, p, t) + res[1] = u[1][1] + res[2] = u[1][2] + res[3] = u[1][3] + res[4] = u[1][4] + res[5] = u[end][2] - 1 + res[6] = u[end][3] - 1 + res[7] = u[end][4] +end +tspan = (0, 700) + +@doc raw""" + flat_earth + +Launch of a satellite into circular orbit from a flat Earth where we assume a uniform gravitational field ``g``. + +Given by + +```math +\frac{dz_1}{dt}=z_3\frac{V_c}{h} +``` +```math +\frac{dz_2}{dt}=z_4\frac{V_c}{h} +``` +```math +\frac{dz_3}{dt}=acc\frac{1}{|V_c|\sqrt{1+z_6^2}} +``` +```math +\frac{dz_4}{dt}=acc\frac{1}{|V_c|\sqrt{1+z_6^2}}-frac{g}{V_c} +``` +```math +\frac{dz_5}{dt}=0 +``` +```math +\frac{dz_6}{dt}=-z_5\frac{V_c}{h} +``` +```math +\frac{dz_7}{dt}=0 +``` + +with boundary condition + +```math +z_1(0)=0 +``` +```math +z_2(0)=0 +``` +```math +z_3(0)=0 +``` +```math +z_4(0)=0 +``` +```math +z_5(1)=h +``` +```math +z_6(1)=V_c +``` +```math +z_7(1)=0 +``` + +# Solution + +No analytical solution + +# References + +[Reference](https://archimede.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=538) +""" +flat_earth = BVProblem(flat_earth_function!, flat_earth_bc!, [0, 0, 0, 0, 0, 0, 0], tspan) + +################### flat_earth_drag ############################ +function flat_earth_drag_function!(du, u, p, t) + fr = 2100000 + h = 180000 + m = 60880 + g_accel = 9.80665 + vc = 1000 * sqrt((398600.4) / (6378.14 + (h / 1000.0))) + beta = 180000 / 840 + eta = 1.225 * 0.5 * 7.069 / 2 + xbardot = u[3] * (vc / h) + ybardot = u[4] * (vc / h) + Vxbardot = (fr / vc * (-u[6] / sqrt(u[6]^2.0 + u[7]^2.0)) - + eta * exp(-u[2] * beta) * u[3] * sqrt(u[3]^2.0 + u[4]^2.0) * vc) / m + Vybardot = (fr / vc * (-u[7] / sqrt(u[6]^2.0 + u[7]^2.0)) - + eta * exp(-u[2] * beta) * u[4] * sqrt(u[3]^2.0 + u[4]^2.0) * vc) / m - + g_accel / vc + if sqrt(u[3]^2.0 + u[4]^2.0) == 0.0 + lambda_2_bar = 0.0 + lambda_3_bar = 0.0 + lambda_4_bar = -u[5] * (vc / h) + else + lambda_2_bar = -(u[6] * u[3] + u[7] * u[4]) * eta * beta * + sqrt(u[3]^2.0 + u[4]^2.0) * exp(-u[2] * beta) * vc / m + lambda_3_bar = eta * exp(-u[2] * beta) * vc * + (u[6] * (2 * u[3]^2.0 + u[4]^2.0) + u[7] * u[3] * u[4]) / + sqrt(u[3]^2.0 + u[4]^2.0) / m + lambda_4_bar = eta * exp(-u[2] * beta) * vc * + (u[7] * (u[3]^2.0 + 2.0 * u[4]^2.0) + u[6] * u[3] * u[4]) / + sqrt(u[3]^2.0 + u[4]^2.0) / m + end + du[1] = xbardot + du[2] = ybardot + du[3] = Vxbardot + du[4] = Vybardot + du[5] = lambda_2_bar + du[6] = lambda_3_bar + du[7] = lambda_4_bar + du[8] = 0 +end + +function flat_earth_drag_bc!(res, u, p, t) + fr = 2100000 + h = 180000 + m = 60880 + g_accel = 9.80665 + vc = 1000 * sqrt((398600.4) / (6378.14 + (h / 1000.0))) + beta = 180000 / 840 + eta = 1.225 * 0.5 * 7.069 / 2 + res[1] = u[1][1] + res[2] = u[1][2] + res[3] = u[1][3] + res[4] = u[1][4] + res[5] = u[end][2] - 1 + res[6] = u[end][3] - 1 + res[7] = u[end][4] + res[8] = (-sqrt(u[6]^2.0 + u[7]^2.0) * fr / m / vc - + (u[6] * u[3]) * eta * exp(-beta) * sqrt(u[3]^2.0) * vc / m - + u[7] * g_accel / vc) * u[8] + 1.0 +end +tspan = (0, 100) +@doc raw""" + flat_earth_drag + +Launch into circular orbit from a flat Earth including athmosferic drag. + +Given by + +```math +\frac{dz_1}{dt}=z_3\frac{V_c}{h} +``` +```math +\frac{dz_2}{dt}=z_4\frac{V_c}{h} +``` +```math +\frac{dz_3}{dt}=\frac{f}{V_c}(-\frac{z_6}{z_6^2+z_7^2}-V_c\eta\exp(-z_2\beta)z_3\sqrt{z_3^3+z_4^2})/m +``` +```math +\frac{dz_4}{dt}=\frac{f}{V_c}(-\frac{z_7}{z_6^2+z_7^2}-V_c\eta\exp(-z_2\beta)z_4\sqrt{z_3^3+z_4^2})/m - g_{accel}/V_c +``` +```math +\frac{dz_5}{dt}=-\eta\beta\exp(-z_2\beta)(z_6z_3+z_7z_4)\sqrt{z_3^3+z_4^2}\frac{V_c}{m} +``` +```math +\frac{dz_6}{dt}=\eta\exp(-z_2\beta)(z_6(2z_3^2+z_4^2)+z_7z_3z_4)V_c/\sqrt{z_3^2+z_4^2}/m +``` +```math +\frac{dz_7}{dt}=\eta\exp(-z_2\beta)(z_7(z_3^2+2z_4^2)+z_6z_3z_4)V_c/\sqrt{z_3^2+z_4^2}/m +``` + +with boundary condition + +```math +z_1(0)=0 +``` +```math +z_2(0)=0 +``` +```math +z_3(0)=0 +``` +```math +z_4(0)=0 +``` +```math +z_5(1)=h +``` +```math +z_6(1)=V_c +``` +```math +z_7(1)=0 +``` + +# Solution + +No analytical solution + +# References + +[Reference](https://archimede.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=544) +""" +flat_earth_drag = BVProblem(flat_earth_drag_function!, + flat_earth_drag_bc!, + [0, 0, 0, 0, 0, 0, 0, 0], + tspan) + +################### measles ############################ +function measles_function!(du, u, p, t) + mu = 0.02 + lambda = 0.0279 + eta = 0.01 + B0 = 1575 + bt = B0 * (1.0 + (cos(2 * π * t))) + + du[1] = mu - (bt * u[1] * u[3]) + du[2] = (bt * u[1] * u[3]) - (u[2] / lambda) + du[3] = (u[2] / lambda) - (u[3] / eta) +end + +function measles_bc!(res, u, p, t) + res[1] = u[1][1] - u[end][1] + res[2] = u[1][2] - u[end][2] + res[3] = u[1][3] - u[end][3] +end +tspan = (0, 1) + +@doc raw""" + measles + +This is an epidemiology model, about the spread of diseases. + +Given by + +```math +\frac{dy_1}{dt}=\mu-\beta(t)y_1y_3 +``` +```math +\frac{dy_2}{dt}=\beta(t)y_1y_3-y_2/\lambda +``` +```math +\frac{dy_3}{dt}=y_2/\lambda-y_3/\eta +``` + +with boundary condition + +```math +y(0)=y(1) +``` + + +# Solution + +No analytical solution + +# References + +[Reference](https://archimede.uniba.it/~bvpsolvers/testsetbvpsolvers/?page_id=555) +""" +measles = BVProblem(measles_function!, measles_bc!, [0, 0, 0], tspan) + # Linear BVP Example Problems export prob_bvp_linear_1, prob_bvp_linear_2, prob_bvp_linear_3, prob_bvp_linear_4, prob_bvp_linear_5, @@ -22,4 +381,6 @@ export prob_bvp_nonlinear_1, prob_bvp_nonlinear_2, prob_bvp_nonlinear_3, prob_bvp_nonlinear_11, prob_bvp_nonlinear_12, prob_bvp_nonlinear_13, prob_bvp_nonlinear_14, prob_bvp_nonlinear_15 +export flat_moon, flat_earth, flat_earth_drag, measles + end # module BVProblemLibrary diff --git a/lib/DAEProblemLibrary/src/DAEProblemLibrary.jl b/lib/DAEProblemLibrary/src/DAEProblemLibrary.jl index c04f071..678f01c 100644 --- a/lib/DAEProblemLibrary/src/DAEProblemLibrary.jl +++ b/lib/DAEProblemLibrary/src/DAEProblemLibrary.jl @@ -34,26 +34,39 @@ Usually solved on ``[0,1e11]`` """ prob_dae_resrob = DAEProblem(f, du0, u0, (0.0, 100000.0)) - -C=[k*1e-6 for k=1:5] -Ub=6; UF=0.026; α=0.99; β=1e-6; R0=1000; R=9000 -Ue(t) = 0.1*sin(200*π*t) -g(x) = 1e-6 * (exp(x/0.026)-1) +C = [k * 1e-6 for k in 1:5] +Ub = 6 +UF = 0.026 +α = 0.99 +β = 1e-6 +R0 = 1000 +R = 9000 +Ue(t) = 0.1 * sin(200 * π * t) +g(x) = 1e-6 * (exp(x / 0.026) - 1) function transamp(out, du, u, p, t) y1, y2, y3, y4, y5, y6, y7, y8 = u - out[1] = -Ue(t)/R0 + y1/R0 + C[1]*du[1] - C[1]*du[2] - out[2] = -Ub/R + y2*2/R - (α-1)*g(y2-y3) - C[1]*du[1] + C[1]*du[2] - out[3] = -g(y2-y3) + y3/R + C[2]*du[3] - out[4] = -Ub/R + y4/R + α*g(y2-y3) + C[3]*du[4] - C[3]*du[5] - out[5] = -Ub/R + y5*2/R - (α-1)*g(y5-y6) - C[3]*du[4] + C[3]*du[5] - out[6] = -g(y5-y6) + y6/R + C[4]*du[6] - out[7] = -Ub/R + y7/R + α*g(y5-y6) + C[5]*du[7] - C[5]*du[8] - out[8] = y8/R - C[5]*du[7] + C[5]*du[8] + out[1] = -Ue(t) / R0 + y1 / R0 + C[1] * du[1] - C[1] * du[2] + out[2] = -Ub / R + y2 * 2 / R - (α - 1) * g(y2 - y3) - C[1] * du[1] + C[1] * du[2] + out[3] = -g(y2 - y3) + y3 / R + C[2] * du[3] + out[4] = -Ub / R + y4 / R + α * g(y2 - y3) + C[3] * du[4] - C[3] * du[5] + out[5] = -Ub / R + y5 * 2 / R - (α - 1) * g(y5 - y6) - C[3] * du[4] + C[3] * du[5] + out[6] = -g(y5 - y6) + y6 / R + C[4] * du[6] + out[7] = -Ub / R + y7 / R + α * g(y5 - y6) + C[5] * du[7] - C[5] * du[8] + out[8] = y8 / R - C[5] * du[7] + C[5] * du[8] end -u0 = [0, Ub/2, Ub/2, Ub, Ub/2, Ub/2, Ub, 0] -du0 = [51.338775, 51.338775, -Ub/(2*(C[2]*R)), -24.9757667, -24.9757667, -Ub/(2*(C[4]*R)), -10.00564453, -10.00564453] +u0 = [0, Ub / 2, Ub / 2, Ub, Ub / 2, Ub / 2, Ub, 0] +du0 = [ + 51.338775, + 51.338775, + -Ub / (2 * (C[2] * R)), + -24.9757667, + -24.9757667, + -Ub / (2 * (C[4] * R)), + -10.00564453, + -10.00564453, +] @doc doc""" The Transistor Amplifier model diff --git a/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl b/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl index 0ee4e26..6b02ecf 100644 --- a/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl +++ b/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl @@ -5,36 +5,36 @@ using DiffEqBase # examples with constant delays export -# DDEs with 1 constant delay - prob_dde_constant_1delay_ip, prob_dde_constant_1delay_oop, - prob_dde_constant_1delay_scalar, - prob_dde_constant_1delay_long_ip, prob_dde_constant_1delay_long_oop, - prob_dde_constant_1delay_long_scalar, -# DDEs with 2 constant delays - prob_dde_constant_2delays_ip, prob_dde_constant_2delays_oop, - prob_dde_constant_2delays_scalar, - prob_dde_constant_2delays_long_ip, prob_dde_constant_2delays_long_oop, - prob_dde_constant_2delays_long_scalar + # DDEs with 1 constant delay + prob_dde_constant_1delay_ip, prob_dde_constant_1delay_oop, + prob_dde_constant_1delay_scalar, + prob_dde_constant_1delay_long_ip, prob_dde_constant_1delay_long_oop, + prob_dde_constant_1delay_long_scalar, + # DDEs with 2 constant delays + prob_dde_constant_2delays_ip, prob_dde_constant_2delays_oop, + prob_dde_constant_2delays_scalar, + prob_dde_constant_2delays_long_ip, prob_dde_constant_2delays_long_oop, + prob_dde_constant_2delays_long_scalar # DDETST problems export -# DDEs with time dependent delays - prob_dde_DDETST_A1, prob_dde_DDETST_A2, -# DDEs with vanishing time dependent delays - prob_dde_DDETST_B1, prob_dde_DDETST_B2, -# DDEs with state dependent delays - prob_dde_DDETST_C1, prob_dde_DDETST_C2, prob_dde_DDETST_C3, prob_dde_DDETST_C4, -# DDEs with vanishing state dependent delays - prob_dde_DDETST_D1, prob_dde_DDETST_D2, -# neutral DDEs with time dependent delays - prob_dde_DDETST_E1, prob_dde_DDETST_E2, -# neutral DDEs with vanishing time dependent delays - prob_dde_DDETST_F1, prob_dde_DDETST_F2, prob_dde_DDETST_F3, prob_dde_DDETST_F4, - prob_dde_DDETST_F5, -# neutral DDEs with state dependent delays - prob_dde_DDETST_G1, prob_dde_DDETST_G2, -# neutral DDEs with vanishing state dependent delays - prob_dde_DDETST_H1, prob_dde_DDETST_H2, prob_dde_DDETST_H3, prob_dde_DDETST_H4 + # DDEs with time dependent delays + prob_dde_DDETST_A1, prob_dde_DDETST_A2, + # DDEs with vanishing time dependent delays + prob_dde_DDETST_B1, prob_dde_DDETST_B2, + # DDEs with state dependent delays + prob_dde_DDETST_C1, prob_dde_DDETST_C2, prob_dde_DDETST_C3, prob_dde_DDETST_C4, + # DDEs with vanishing state dependent delays + prob_dde_DDETST_D1, prob_dde_DDETST_D2, + # neutral DDEs with time dependent delays + prob_dde_DDETST_E1, prob_dde_DDETST_E2, + # neutral DDEs with vanishing time dependent delays + prob_dde_DDETST_F1, prob_dde_DDETST_F2, prob_dde_DDETST_F3, prob_dde_DDETST_F4, + prob_dde_DDETST_F5, + # neutral DDEs with state dependent delays + prob_dde_DDETST_G1, prob_dde_DDETST_G2, + # neutral DDEs with vanishing state dependent delays + prob_dde_DDETST_H1, prob_dde_DDETST_H2, prob_dde_DDETST_H3, prob_dde_DDETST_H4 # RADAR5 problems export prob_dde_RADAR5_oregonator, prob_dde_RADAR5_robertson, prob_dde_RADAR5_waltman diff --git a/lib/DDEProblemLibrary/src/constant_delays.jl b/lib/DDEProblemLibrary/src/constant_delays.jl index d3918a6..7ac4914 100644 --- a/lib/DDEProblemLibrary/src/constant_delays.jl +++ b/lib/DDEProblemLibrary/src/constant_delays.jl @@ -16,7 +16,7 @@ but not of `u`). """ function remake_dde_constant_u0_tType(prob::DDEProblem, u₀, tType) remake(prob; u0 = u₀, tspan = tType.(prob.tspan), p = eltype(u₀), - constant_lags = tType.(prob.constant_lags)) + constant_lags = tType.(prob.constant_lags)) end # History functions @@ -72,10 +72,10 @@ function f_dde_constant_1delay_ip!(du, u, h, p, t) end function fanalytic_dde_constant_1delay(u₀, - ::Union{typeof(h_dde_constant_ip), - typeof(h_dde_constant_oop), - typeof(h_dde_constant_scalar)}, - p, t) + ::Union{typeof(h_dde_constant_ip), + typeof(h_dde_constant_oop), + typeof(h_dde_constant_scalar)}, + p, t) z = t * inv(oneunit(t)) 0 ≤ z ≤ 10 || error("analytical solution is only implemented for t ∈ [0, 10]") @@ -96,14 +96,14 @@ function fanalytic_dde_constant_1delay(u₀, c = @evalpoly(z, 13201//120, -13081//120, 521//12, -107//12, 1, -7//120, 1//720) elseif z < 8 c = @evalpoly(z, 39371//144, -39227//144, 27227//240, -3685//144, 487//144, - -21//80, 1//90, -1//5040) + -21//80, 1//90, -1//5040) elseif z < 9 c = @evalpoly(z, 1158379//1680, -1156699//1680, 212753//720, -51193//720, - 1511//144, -701//720, 1//18, -1//560, 1//40320) + 1511//144, -701//720, 1//18, -1//560, 1//40320) else c = @evalpoly(z, 23615939//13440, -23602499//13440, 7761511//10080, - -279533//1440, 89269//2880, -1873//576, 323//1440, -11//1120, - 1//4032, -1//362880) + -279533//1440, 89269//2880, -1873//576, 323//1440, -11//1120, + 1//4032, -1//362880) end c .* u₀ @@ -111,10 +111,10 @@ function fanalytic_dde_constant_1delay(u₀, end const prob_dde_constant_1delay_ip = DDEProblem(DDEFunction(f_dde_constant_1delay_ip!; - analytic = fanalytic_dde_constant_1delay), - [1.0], h_dde_constant_ip, (0.0, 10.0), - typeof(1.0); - constant_lags = [1]) + analytic = fanalytic_dde_constant_1delay), + [1.0], h_dde_constant_ip, (0.0, 10.0), + typeof(1.0); + constant_lags = [1]) ### Out-of-place function @@ -132,10 +132,10 @@ function f_dde_constant_1delay_oop(u, h, p, t) end const prob_dde_constant_1delay_oop = DDEProblem(DDEFunction(f_dde_constant_1delay_oop; - analytic = fanalytic_dde_constant_1delay), - [1.0], h_dde_constant_oop, (0.0, 10.0), - typeof(1.0); - constant_lags = [1]) + analytic = fanalytic_dde_constant_1delay), + [1.0], h_dde_constant_oop, (0.0, 10.0), + typeof(1.0); + constant_lags = [1]) ### Scalar function @@ -153,10 +153,10 @@ function f_dde_constant_1delay_scalar(u, h, p, t) end const prob_dde_constant_1delay_scalar = DDEProblem(DDEFunction(f_dde_constant_1delay_scalar; - analytic = fanalytic_dde_constant_1delay), - 1.0, h_dde_constant_scalar, (0.0, 10.0), - typeof(1.0); - constant_lags = [1]) + analytic = fanalytic_dde_constant_1delay), + 1.0, h_dde_constant_scalar, (0.0, 10.0), + typeof(1.0); + constant_lags = [1]) ## Long time span @@ -184,9 +184,9 @@ function f_dde_constant_1delay_long_ip!(du, u, h, p, t) end const prob_dde_constant_1delay_long_ip = DDEProblem(f_dde_constant_1delay_long_ip!, [1.0], - h_dde_constant_ip, (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 5]) + h_dde_constant_ip, (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 5]) ### Out-of-place function @@ -204,9 +204,9 @@ function f_dde_constant_1delay_long_oop(u, h, p, t) end const prob_dde_constant_1delay_long_oop = DDEProblem(f_dde_constant_1delay_long_oop, [1.0], - h_dde_constant_oop, (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 5]) + h_dde_constant_oop, (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 5]) ### Scalar function @@ -224,10 +224,10 @@ function f_dde_constant_1delay_long_scalar(u, h, p, t) end const prob_dde_constant_1delay_long_scalar = DDEProblem(f_dde_constant_1delay_long_scalar, - 1.0, h_dde_constant_scalar, - (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 5]) + 1.0, h_dde_constant_scalar, + (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 5]) # Two constant delays @@ -261,10 +261,10 @@ function f_dde_constant_2delays_ip!(du, u, h, p, t) end function fanalytic_dde_constant_2delays(u₀, - ::Union{typeof(h_dde_constant_ip), - typeof(h_dde_constant_oop), - typeof(h_dde_constant_scalar)}, - p, t) + ::Union{typeof(h_dde_constant_ip), + typeof(h_dde_constant_oop), + typeof(h_dde_constant_scalar)}, + p, t) z = t * inv(oneunit(t)) 0 ≤ z ≤ 1 || error("analytical solution is only implemented for t ∈ [0, 1]") @@ -298,10 +298,10 @@ function fanalytic_dde_constant_2delays(u₀, end const prob_dde_constant_2delays_ip = DDEProblem(DDEFunction(f_dde_constant_2delays_ip!; - analytic = fanalytic_dde_constant_2delays), - [1.0], h_dde_constant_ip, (0.0, 1.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + analytic = fanalytic_dde_constant_2delays), + [1.0], h_dde_constant_ip, (0.0, 1.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) ### Out-of-place function @@ -319,10 +319,10 @@ function f_dde_constant_2delays_oop(u, h, p, t) end const prob_dde_constant_2delays_oop = DDEProblem(DDEFunction(f_dde_constant_2delays_oop; - analytic = fanalytic_dde_constant_2delays), - [1.0], h_dde_constant_oop, (0.0, 1.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + analytic = fanalytic_dde_constant_2delays), + [1.0], h_dde_constant_oop, (0.0, 1.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) ### Scalar function @@ -340,10 +340,10 @@ function f_dde_constant_2delays_scalar(u, h, p, t) end const prob_dde_constant_2delays_scalar = DDEProblem(DDEFunction(f_dde_constant_2delays_scalar; - analytic = fanalytic_dde_constant_2delays), - 1.0, h_dde_constant_scalar, (0.0, 1.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + analytic = fanalytic_dde_constant_2delays), + 1.0, h_dde_constant_scalar, (0.0, 1.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) ## Long time span @@ -370,9 +370,9 @@ function f_dde_constant_2delays_long_ip!(du, u, h, p, t) end const prob_dde_constant_2delays_long_ip = DDEProblem(f_dde_constant_2delays_long_ip!, [1.0], - h_dde_constant_ip, (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + h_dde_constant_ip, (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) ### Out-of-place function @@ -390,10 +390,10 @@ function f_dde_constant_2delays_long_oop(u, h, p, t) end const prob_dde_constant_2delays_long_oop = DDEProblem(f_dde_constant_2delays_long_oop, - [1.0], h_dde_constant_oop, - (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + [1.0], h_dde_constant_oop, + (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) #### Scalar function @@ -411,7 +411,7 @@ function f_dde_constant_2delays_long_scalar(u, h, p, t) end const prob_dde_constant_2delays_long_scalar = DDEProblem(f_dde_constant_2delays_long_scalar, - 1.0, h_dde_constant_scalar, - (0.0, 100.0), - typeof(1.0); - constant_lags = [1 / 3, 1 / 5]) + 1.0, h_dde_constant_scalar, + (0.0, 100.0), + typeof(1.0); + constant_lags = [1 / 3, 1 / 5]) diff --git a/lib/DDEProblemLibrary/src/ddetst.jl b/lib/DDEProblemLibrary/src/ddetst.jl index 63b32a8..d829d07 100644 --- a/lib/DDEProblemLibrary/src/ddetst.jl +++ b/lib/DDEProblemLibrary/src/ddetst.jl @@ -36,7 +36,7 @@ function h_dde_DDETST_A1(p, t) end const prob_dde_DDETST_A1 = DDEProblem(f_dde_DDETST_A1, h_dde_DDETST_A1, (0.0, 500.0); - constant_lags = [14]) + constant_lags = [14]) # Problem A2 @doc raw""" @@ -93,7 +93,7 @@ function h_dde_DDETST_A2(p, t; idxs::Union{Nothing, Int} = nothing) end const prob_dde_DDETST_A2 = DDEProblem(f_dde_DDETST_A2!, h_dde_DDETST_A2, (0.0, 100.0); - constant_lags = [20]) + constant_lags = [20]) # Problem B2 @doc raw""" @@ -140,9 +140,9 @@ function fanalytic_dde_DDETST_B1(u₀, ::typeof(h_dde_DDETST_B1), p, t) end const prob_dde_DDETST_B1 = DDEProblem(DDEFunction(f_dde_DDETST_B1; - analytic = fanalytic_dde_DDETST_B1), - h_dde_DDETST_B1, (0.1, 10.0); - dependent_lags = ((u, p, t) -> t - exp(1 - 1 / t),)) + analytic = fanalytic_dde_DDETST_B1), + h_dde_DDETST_B1, (0.1, 10.0); + dependent_lags = ((u, p, t) -> t - exp(1 - 1 / t),)) # Problem B2 @doc raw""" @@ -204,9 +204,9 @@ function fanalytic_dde_DDETST_B2(u₀, ::typeof(h_dde_DDETST_B2), p, t) end const prob_dde_DDETST_B2 = DDEProblem(DDEFunction(f_dde_DDETST_B2; - analytic = fanalytic_dde_DDETST_B2), - h_dde_DDETST_B2, (0.0, 2 * log(66)); - dependent_lags = ((u, p, t) -> t / 2,)) + analytic = fanalytic_dde_DDETST_B2), + h_dde_DDETST_B2, (0.0, 2 * log(66)); + dependent_lags = ((u, p, t) -> t / 2,)) # Problem C1 @doc raw""" @@ -236,7 +236,7 @@ function h_dde_DDETST_C1(p, t) end const prob_dde_DDETST_C1 = DDEProblem(f_dde_DDETST_C1, h_dde_DDETST_C1, (0.0, 30.0); - dependent_lags = ((u, p, t) -> 1 + abs(u),)) + dependent_lags = ((u, p, t) -> 1 + abs(u),)) # Problem C2 @doc raw""" @@ -294,7 +294,7 @@ function h_dde_DDETST_C2(p, t; idxs::Union{Nothing, Int} = nothing) end const prob_dde_DDETST_C2 = DDEProblem(f_dde_DDETST_C2!, h_dde_DDETST_C2, (0.0, 30.0); - dependent_lags = ((u, p, t) -> u[2],)) + dependent_lags = ((u, p, t) -> u[2],)) # Problem C3 @doc raw""" @@ -357,7 +357,7 @@ const prob_dde_DDETST_C3 = let s₀ = 0.0031, T₁ = 6, γ = 0.001, Q = 0.0275, end DDEProblem(f_dde_DDETST_C3!, h_dde_DDETST_C3, (0.0, 300.0); - constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],)) + constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],)) end # Problem C4 @@ -415,7 +415,7 @@ const prob_dde_DDETST_C4 = let s₀ = 0.00372, T₁ = 3, γ = 0.01, Q = 0.00178, end DDEProblem(f_dde_DDETST_C4!, [3.5, 10.0, 50.0], h_dde_DDETST_C4, (0.0, 100.0); - constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],)) + constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],)) end # Problem D1 @@ -491,9 +491,9 @@ function fanalytic_dde_DDETST_D1(u₀, ::typeof(h_dde_DDETST_D1), p, t) end const prob_dde_DDETST_D1 = DDEProblem(DDEFunction(f_dde_DDETST_D1!; - analytic = fanalytic_dde_DDETST_D1), - h_dde_DDETST_D1, (0.1, 5.0); - dependent_lags = ((u, p, t) -> t - exp(1 - u[2]),)) + analytic = fanalytic_dde_DDETST_D1), + h_dde_DDETST_D1, (0.1, 5.0); + dependent_lags = ((u, p, t) -> t - exp(1 - u[2]),)) # Problem D2 @doc raw""" @@ -561,7 +561,7 @@ const prob_dde_DDETST_D2 = let r₁ = 0.02, r₂ = 0.005, α = 3, δ = 0.01 end DDEProblem(f_dde_DDETST_D2!, h_dde_DDETST_D2, (0.0, 40.0); - dependent_lags = ((u, p, t) -> u[4],)) + dependent_lags = ((u, p, t) -> u[4],)) end # Problem E1 @@ -598,7 +598,7 @@ const prob_dde_DDETST_E1 = let r = π / sqrt(3) + 1 / 20, c = sqrt(3) / (2 * π) end DDEProblem(f_dde_DDETST_E1, h_dde_DDETST_E1, (0.0, 40.0); - constant_lags = [1], neutral = true) + constant_lags = [1], neutral = true) end # Problem E2 @@ -656,7 +656,7 @@ const prob_dde_DDETST_E2 = let α = 0.1, ρ = 2.9, τ = 0.42 end global function h_dde_DDETST_E2(p, t, ::Type{Val{1}}; - idxs::Union{Nothing, Int} = nothing) + idxs::Union{Nothing, Int} = nothing) t ≤ 0 || error("history function is only implemented for t ≤ 0") if idxs === nothing @@ -669,7 +669,7 @@ const prob_dde_DDETST_E2 = let α = 0.1, ρ = 2.9, τ = 0.42 end DDEProblem(f_dde_DDETST_E2!, h_dde_DDETST_E2, (0.0, 2.0); - constant_lags = [τ], neutral = true) + constant_lags = [τ], neutral = true) end # Problem F1 @@ -730,10 +730,10 @@ function fanalytic_dde_DDETST_F1(u₀, ::typeof(h_dde_DDETST_F1), p, t) end const prob_dde_DDETST_F1 = DDEProblem(DDEFunction(f_dde_DDETST_F1; - analytic = fanalytic_dde_DDETST_F1), - h_dde_DDETST_F1, (0.0, 0.1); - dependent_lags = ((u, p, t) -> t / 2,), - neutral = true) + analytic = fanalytic_dde_DDETST_F1), + h_dde_DDETST_F1, (0.0, 0.1); + dependent_lags = ((u, p, t) -> t / 2,), + neutral = true) # Problem F2 @doc raw""" @@ -820,10 +820,10 @@ function fanalytic_dde_DDETST_F2(u₀, ::typeof(h_dde_DDETST_F2), p, t) end const prob_dde_DDETST_F2 = DDEProblem(DDEFunction(f_dde_DDETST_F2; - analytic = fanalytic_dde_DDETST_F2), - h_dde_DDETST_F2, (0.25, 0.499); - dependent_lags = ((u, p, t) -> 1 / 2 - t,), - neutral = true) + analytic = fanalytic_dde_DDETST_F2), + h_dde_DDETST_F2, (0.25, 0.499); + dependent_lags = ((u, p, t) -> 1 / 2 - t,), + neutral = true) # Problem F3 @doc raw""" @@ -876,11 +876,11 @@ function fanalytic_dde_DDETST_F345(u₀, ::typeof(h_dde_DDETST_F345), p, t) end const prob_dde_DDETST_F3 = DDEProblem(DDEFunction(f_dde_DDETST_F3; - analytic = fanalytic_dde_DDETST_F345), - h_dde_DDETST_F345, (0.0, 10.0); - dependent_lags = ((u, p, t) -> 0.5 * t * - (1 + cos(2 * π * t)),), - neutral = true) + analytic = fanalytic_dde_DDETST_F345), + h_dde_DDETST_F345, (0.0, 10.0); + dependent_lags = ((u, p, t) -> 0.5 * t * + (1 + cos(2 * π * t)),), + neutral = true) # Problem F4 """ @@ -899,8 +899,8 @@ let L₃ = 0.4 end const prob_dde_DDETST_F4 = remake(prob_dde_DDETST_F3; - f = DDEFunction(f_dde_DDETST_F4; - analytic = fanalytic_dde_DDETST_F345)) + f = DDEFunction(f_dde_DDETST_F4; + analytic = fanalytic_dde_DDETST_F345)) # Problem F5 """ @@ -919,8 +919,8 @@ let L₃ = 0.6 end const prob_dde_DDETST_F5 = remake(prob_dde_DDETST_F3; - f = DDEFunction(f_dde_DDETST_F5; - analytic = fanalytic_dde_DDETST_F345)) + f = DDEFunction(f_dde_DDETST_F5; + analytic = fanalytic_dde_DDETST_F345)) # Problem G1 @doc raw""" @@ -972,10 +972,10 @@ function fanalytic_dde_DDETST_G1(u₀, ::typeof(h_dde_DDETST_G1), p, t) end const prob_dde_DDETST_G1 = DDEProblem(DDEFunction(f_dde_DDETST_G1; - analytic = fanalytic_dde_DDETST_G1), - h_dde_DDETST_G1, (0.0, 1.0); - dependent_lags = ((u, p, t) -> u^2 / 4,), - neutral = true) + analytic = fanalytic_dde_DDETST_G1), + h_dde_DDETST_G1, (0.0, 1.0); + dependent_lags = ((u, p, t) -> u^2 / 4,), + neutral = true) # Problem G2 @doc raw""" @@ -1025,10 +1025,10 @@ function fanalytic_dde_DDETST_G2(u₀, ::typeof(h_dde_DDETST_G2), p, t) end const prob_dde_DDETST_G2 = DDEProblem(DDEFunction(f_dde_DDETST_G2; - analytic = fanalytic_dde_DDETST_G2), - h_dde_DDETST_G2, (0.0, 1.0); - dependent_lags = ((u, p, t) -> t + 2 - u,), - neutral = true) + analytic = fanalytic_dde_DDETST_G2), + h_dde_DDETST_G2, (0.0, 1.0); + dependent_lags = ((u, p, t) -> t + 2 - u,), + neutral = true) # Problem H1 @doc raw""" @@ -1084,10 +1084,10 @@ function fanalytic_dde_DDETST_H1(u₀, ::typeof(h_dde_DDETST_H1), p, t) end const prob_dde_DDETST_H1 = DDEProblem(DDEFunction(f_dde_DDETST_H1; - analytic = fanalytic_dde_DDETST_H1), - h_dde_DDETST_H1, (0.0, 0.225 * π); - dependent_lags = ((u, p, t) -> t / (1 + u^2),), - neutral = true) + analytic = fanalytic_dde_DDETST_H1), + h_dde_DDETST_H1, (0.0, 0.225 * π); + dependent_lags = ((u, p, t) -> t / (1 + u^2),), + neutral = true) # Problem H2 @doc raw""" @@ -1149,10 +1149,10 @@ function fanalytic_dde_DDETST_H234(u₀, ::typeof(h_dde_DDETST_H234), p, t) end const prob_dde_DDETST_H2 = DDEProblem(DDEFunction(f_dde_DDETST_H2; - analytic = fanalytic_dde_DDETST_H234), - h_dde_DDETST_H234, (0.0, π); - dependent_lags = ((u, p, t) -> t * (1 - u^2),), - neutral = true) + analytic = fanalytic_dde_DDETST_H234), + h_dde_DDETST_H234, (0.0, π); + dependent_lags = ((u, p, t) -> t * (1 - u^2),), + neutral = true) # Problem H3 """ @@ -1180,8 +1180,8 @@ let L₃ = 0.3 end const prob_dde_DDETST_H3 = remake(prob_dde_DDETST_H2; - f = DDEFunction(f_dde_DDETST_H3; - analytic = fanalytic_dde_DDETST_H234)) + f = DDEFunction(f_dde_DDETST_H3; + analytic = fanalytic_dde_DDETST_H234)) # Problem H4 """ @@ -1209,5 +1209,5 @@ let L₃ = 0.5 end const prob_dde_DDETST_H4 = remake(prob_dde_DDETST_H2; - f = DDEFunction(f_dde_DDETST_H4; - analytic = fanalytic_dde_DDETST_H234)) + f = DDEFunction(f_dde_DDETST_H4; + analytic = fanalytic_dde_DDETST_H234)) diff --git a/lib/DDEProblemLibrary/src/radar5.jl b/lib/DDEProblemLibrary/src/radar5.jl index 15dc421..7f7a04c 100644 --- a/lib/DDEProblemLibrary/src/radar5.jl +++ b/lib/DDEProblemLibrary/src/radar5.jl @@ -64,7 +64,7 @@ const prob_dde_RADAR5_oregonator = let k₁ = 1.34, k₂ = 1.6e9, k₃ = 8_000, end DDEProblem(f_dde_RADAR5_oregonator!, h_dde_RADAR5_oregonator, (0.0, 100.5); - constant_lags = [τ]) + constant_lags = [τ]) end # ROBERTSON example @@ -122,7 +122,7 @@ const prob_dde_RADAR5_robertson = let a = 0.04, b = 10_000, c = 3e7, τ = 0.01 end DDEProblem(f_dde_RADAR5_robertson!, h_dde_RADAR5_robertson, (0.0, 1e10); - constant_lags = [τ]) + constant_lags = [τ]) end # WALTMAN example @@ -224,12 +224,12 @@ function h_dde_RADAR5_waltman(p, t) end const prob_dde_RADAR5_waltman = DDEProblem(f_dde_RADAR5_waltman!, - (p, t) -> h_dde_RADAR5_waltman(p, 0.0), - h_dde_RADAR5_waltman, (0.0, 300.0), - (ϕ₀ = 0.75e-4, t₀ = 32, t₁ = 119); - dependent_lags = ((u, p, t) -> t - u[5], - (u, p, t) -> t - u[6]), - order_discontinuity_t0 = 1) + (p, t) -> h_dde_RADAR5_waltman(p, 0.0), + h_dde_RADAR5_waltman, (0.0, 300.0), + (ϕ₀ = 0.75e-4, t₀ = 32, t₁ = 119); + dependent_lags = ((u, p, t) -> t - u[5], + (u, p, t) -> t - u[6]), + order_discontinuity_t0 = 1) const prob_dde_RADAR5_waltman_1 = prob_dde_RADAR5_waltman @doc raw""" @@ -244,7 +244,7 @@ Waltman, P. (1978). A threshold model of antigen-stimulated antibody production, Immunology (8), pp. 437-453. """ const prob_dde_RADAR5_waltman_2 = remake(prob_dde_RADAR5_waltman; - p = (ϕ₀ = 0.5e-4, t₀ = 32, t₁ = 111)) + p = (ϕ₀ = 0.5e-4, t₀ = 32, t₁ = 111)) @doc raw""" prob_dde_RADAR5_waltman_3 @@ -258,7 +258,7 @@ Waltman, P. (1978). A threshold model of antigen-stimulated antibody production, Immunology (8), pp. 437-453. """ const prob_dde_RADAR5_waltman_3 = remake(prob_dde_RADAR5_waltman; - p = (ϕ₀ = 1e-5, t₀ = 33, t₁ = 145)) + p = (ϕ₀ = 1e-5, t₀ = 33, t₁ = 145)) @doc raw""" prob_dde_RADAR5_waltman_4 @@ -272,7 +272,7 @@ Waltman, P. (1978). A threshold model of antigen-stimulated antibody production, Immunology (8), pp. 437-453. """ const prob_dde_RADAR5_waltman_4 = remake(prob_dde_RADAR5_waltman; - p = (ϕ₀ = 0.75e-5, t₀ = 34, t₁ = 163)) + p = (ϕ₀ = 0.75e-5, t₀ = 34, t₁ = 163)) @doc raw""" prob_dde_RADAR5_waltman_5 @@ -286,4 +286,4 @@ Waltman, P. (1978). A threshold model of antigen-stimulated antibody production, Immunology (8), pp. 437-453. """ const prob_dde_RADAR5_waltman_5 = remake(prob_dde_RADAR5_waltman; - p = (ϕ₀ = 0.5e-5, t₀ = 35, t₁ = 197)) + p = (ϕ₀ = 0.5e-5, t₀ = 35, t₁ = 197)) diff --git a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl index 088e458..2f82a08 100644 --- a/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl +++ b/lib/JumpProblemLibrary/src/JumpProblemLibrary.jl @@ -9,12 +9,12 @@ RuntimeGeneratedFunctions.init(@__MODULE__) # Jump Example Problems export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs, -# examples mixing mass action and constant rate jumps - prob_jump_osc_mixed_jumptypes, -# examples used in published benchmarks / comparisions - prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor, -# examples approximating diffusion by continuous time random walks - prob_jump_diffnetwork + # examples mixing mass action and constant rate jumps + prob_jump_osc_mixed_jumptypes, + # examples used in published benchmarks / comparisions + prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor, + # examples approximating diffusion by continuous time random walks + prob_jump_diffnetwork """ General structure to hold JumpProblem info. Needed since @@ -107,16 +107,16 @@ prob = DiscreteProblem(oscil_rs, u0, (0.0, tf), eval_module = @__MODULE__) prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, prob, nothing) specs_sym_to_name = Dict(:S1 => "R(a,l)", - :S2 => "L(r)", - :S3 => "A(Y~U,r)", - :S4 => "L(r!1).R(a,l!1)", - :S5 => "A(Y~U,r!1).R(a!1,l)", - :S6 => "A(Y~U,r!1).L(r!2).R(a!1,l!2)", - :S7 => "A(Y~P,r!1).L(r!2).R(a!1,l!2)", - :S8 => "A(Y~P,r!1).R(a!1,l)", - :S9 => "A(Y~P,r)") + :S2 => "L(r)", + :S3 => "A(Y~U,r)", + :S4 => "L(r!1).R(a,l!1)", + :S5 => "A(Y~U,r!1).R(a!1,l)", + :S6 => "A(Y~U,r!1).L(r!2).R(a!1,l!2)", + :S7 => "A(Y~P,r!1).L(r!2).R(a!1,l!2)", + :S8 => "A(Y~P,r!1).R(a!1,l)", + :S9 => "A(Y~P,r)") rates_sym_to_idx = Dict(:R0 => 1, :L0 => 2, :A0 => 3, :kon => 4, :koff => 5, - :kAon => 6, :kAoff => 7, :kAp => 8, :kAdp => 9) + :kAon => 6, :kAoff => 7, :kAp => 8, :kAdp => 9) params = [5360, 1160, 5360, 0.01, 0.1, 0.01, 0.1, 0.01, 0.1] rs = @reaction_network begin kon, S1 + S2 --> S4 @@ -154,9 +154,9 @@ prob = DiscreteProblem(rs, u0, (0.0, tf), rates, eval_module = @__MODULE__) Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml """ prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob, - Dict("specs_to_sym_name" => specs_sym_to_name, - "rates_sym_to_idx" => rates_sym_to_idx, - "params" => params)) + Dict("specs_to_sym_name" => specs_sym_to_name, + "rates_sym_to_idx" => rates_sym_to_idx, + "params" => params)) # generate the network N = 10 # number of genes @@ -170,9 +170,9 @@ function construct_genenetwork(N) addspecies!(genenetwork, M[2 * i - i]) addspecies!(genenetwork, P[2 * i - i]) addreaction!(genenetwork, - Reaction(10.0, [G[2 * i - i]], [G[2 * i - i], M[2 * i - i]])) + Reaction(10.0, [G[2 * i - i]], [G[2 * i - i], M[2 * i - i]])) addreaction!(genenetwork, - Reaction(10.0, [M[2 * i - i]], [M[2 * i - i], P[2 * i - i]])) + Reaction(10.0, [M[2 * i - i]], [M[2 * i - i], P[2 * i - i]])) addreaction!(genenetwork, Reaction(1.0, [M[2 * i - i]], nothing)) addreaction!(genenetwork, Reaction(1.0, [P[2 * i - i]], nothing)) # genenetwork *= "\t 10.0, G$(2*i-1) --> G$(2*i-1) + M$(2*i-1)\n" @@ -194,7 +194,7 @@ function construct_genenetwork(N) addspecies!(genenetwork, G_ind[2 * i]) addreaction!(genenetwork, - Reaction(0.0001, [G[2 * i], P[2 * i - i]], [G_ind[2 * i]])) + Reaction(0.0001, [G[2 * i], P[2 * i - i]], [G_ind[2 * i]])) addreaction!(genenetwork, Reaction(100.0, [G_ind[2 * i]], [G_ind[2 * i], M[2 * i]])) # genenetwork *= "\t 0.0001, G$(2*i) + P$(2*i-1) --> G$(2*i)_ind \n" # genenetwork *= "\t 100., G$(2*i)_ind --> G$(2*i)_ind + M$(2*i)\n" @@ -238,7 +238,7 @@ prob = DiscreteProblem(rn, u0, (0.0, tf), rnpar, eval_module = @__MODULE__) Springer (2017). """ prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob, - Dict("specs_names" => varlabels)) + Dict("specs_names" => varlabels)) # diffusion model function getDiffNetwork(N) @@ -267,6 +267,6 @@ tf = 10.0 u0 is a similar function that returns the initial condition vector. """ prob_jump_diffnetwork = JumpProblemNetwork(getDiffNetwork, params, tf, getDiffu0, nothing, - nothing) + nothing) end # module diff --git a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl index c454769..6d5d8b5 100644 --- a/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl +++ b/lib/ODEProblemLibrary/src/ODEProblemLibrary.jl @@ -17,19 +17,19 @@ Random.seed!(100) #ODE Example Problems export prob_ode_linear, prob_ode_bigfloatlinear, prob_ode_2Dlinear, - prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody, - prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff, - prob_ode_lotkavolterra, prob_ode_fitzhughnagumo, - prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades, - prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego, - prob_ode_hires, prob_ode_pollution, prob_ode_filament, - prob_ode_nonlinchem, - SolverDiffEq, filament_prob, - thomas, lorenz, aizawa, dadras, chen, rossler, rabinovich_fabrikant, sprott, - hindmarsh_rose, - prob_ode_thomas, prob_ode_lorenz, prob_ode_aizawa, prob_ode_dadras, - prob_ode_chen, prob_ode_rossler, prob_ode_rabinovich_fabrikant, prob_ode_sprott, - prob_ode_hindmarsh_rose + prob_ode_large2Dlinear, prob_ode_bigfloat2Dlinear, prob_ode_rigidbody, + prob_ode_2Dlinear_notinplace, prob_ode_vanderpol, prob_ode_vanderpol_stiff, + prob_ode_lotkavolterra, prob_ode_fitzhughnagumo, + prob_ode_rober, prob_ode_threebody, prob_ode_mm_linear, prob_ode_pleiades, + prob_ode_brusselator_1d, prob_ode_brusselator_2d, prob_ode_orego, + prob_ode_hires, prob_ode_pollution, prob_ode_filament, + prob_ode_nonlinchem, + SolverDiffEq, filament_prob, + thomas, lorenz, aizawa, dadras, chen, rossler, rabinovich_fabrikant, sprott, + hindmarsh_rose, + prob_ode_thomas, prob_ode_lorenz, prob_ode_aizawa, prob_ode_dadras, + prob_ode_chen, prob_ode_rossler, prob_ode_rabinovich_fabrikant, prob_ode_sprott, + prob_ode_hindmarsh_rose include("ode_linear_prob.jl") include("ode_simple_nonlinear_prob.jl") diff --git a/lib/ODEProblemLibrary/src/brusselator_prob.jl b/lib/ODEProblemLibrary/src/brusselator_prob.jl index a0d3ca2..cbaca33 100644 --- a/lib/ODEProblemLibrary/src/brusselator_prob.jl +++ b/lib/ODEProblemLibrary/src/brusselator_prob.jl @@ -1,6 +1,6 @@ function brusselator_f(x, y, t) ifelse((((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) && - (t >= 1.1), 5.0, 0.0) + (t >= 1.1), 5.0, 0.0) end function limit(a, N) if a == N + 1 @@ -90,11 +90,11 @@ u(x,y+1,t) = u(x,y,t) From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 152 """ prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, - init_brusselator_2d(xyd_brusselator), - (0.0, 11.5), - (3.4, 1.0, 10.0, - xyd_brusselator, step(xyd_brusselator), - length(xyd_brusselator))) + init_brusselator_2d(xyd_brusselator), + (0.0, 11.5), + (3.4, 1.0, 10.0, + xyd_brusselator, step(xyd_brusselator), + length(xyd_brusselator))) const N_brusselator_1d = 40 @@ -152,6 +152,6 @@ v(0,t) = v(1,t) = 3 From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 6 """ prob_ode_brusselator_1d = ODEProblem(brusselator_1d_loop, - init_brusselator_1d(N_brusselator_1d), - (0.0, 10.0), - (1.0, 3.0, 1 / 41, zeros(N_brusselator_1d))) + init_brusselator_1d(N_brusselator_1d), + (0.0, 10.0), + (1.0, 3.0, 1 / 41, zeros(N_brusselator_1d))) diff --git a/lib/ODEProblemLibrary/src/filament_prob.jl b/lib/ODEProblemLibrary/src/filament_prob.jl index ae7b7b2..a6239df 100644 --- a/lib/ODEProblemLibrary/src/filament_prob.jl +++ b/lib/ODEProblemLibrary/src/filament_prob.jl @@ -10,10 +10,10 @@ struct FerromagneticContinuous <: AbstractMagneticForce end mutable struct FilamentCache{ - MagneticForce <: AbstractMagneticForce, - InextensibilityCache <: AbstractInextensibilityCache, - SolverCache <: AbstractSolverCache - } <: AbstractFilamentCache + MagneticForce <: AbstractMagneticForce, + InextensibilityCache <: AbstractInextensibilityCache, + SolverCache <: AbstractSolverCache, +} <: AbstractFilamentCache N::Int μ::T Cm::T @@ -36,7 +36,8 @@ struct NoHydroProjectionCache <: AbstractInextensibilityCache new(zeros(N, 3 * (N + 1)), # J zeros(3 * (N + 1), 3 * (N + 1)), # P zeros(N, N), # J_JT - LinearAlgebra.LDLt{T, SymTridiagonal{T}}(SymTridiagonal(zeros(N), zeros(N - 1))), + LinearAlgebra.LDLt{T, SymTridiagonal{T}}(SymTridiagonal(zeros(N), + zeros(N - 1))), zeros(N, 3 * (N + 1))) end end @@ -51,22 +52,22 @@ function FilamentCache(N = 20; Cm = 32, ω = 200, Solver = SolverDiffEq) SolverCache = DiffEqSolverCache tmp = zeros(3 * (N + 1)) FilamentCache{FerromagneticContinuous, InextensibilityCache, SolverCache}(N, N + 1, Cm, - view(tmp, - 1:3:(3 * (N + 1))), - view(tmp, - 2:3:(3 * (N + 1))), - view(tmp, - 3:3:(3 * (N + 1))), - zeros(3 * - (N + 1), - 3 * - (N + 1)), # A - InextensibilityCache(N), # P - FerromagneticContinuous(ω, - zeros(3 * - (N + - 1))), - SolverCache(N)) + view(tmp, + 1:3:(3 * (N + 1))), + view(tmp, + 2:3:(3 * (N + 1))), + view(tmp, + 3:3:(3 * (N + 1))), + zeros(3 * + (N + 1), + 3 * + (N + 1)), # A + InextensibilityCache(N), # P + FerromagneticContinuous(ω, + zeros(3 * + (N + + 1))), + SolverCache(N)) end function stiffness_matrix!(f::AbstractFilamentCache) N, μ, A = f.N, f.μ, f.A @@ -204,7 +205,7 @@ function subtract_from_identity!(A) end function LDLt_inplace!(L::LinearAlgebra.LDLt{T, SymTridiagonal{T}}, - A::Matrix{T}) where {T <: Real} + A::Matrix{T}) where {T <: Real} n = size(A, 1) dv, ev = L.data.dv, L.data.ev @inbounds for (i, d) in enumerate(diagind(A)) diff --git a/lib/ODEProblemLibrary/src/nonlinchem.jl b/lib/ODEProblemLibrary/src/nonlinchem.jl index 54415a7..003d0e0 100644 --- a/lib/ODEProblemLibrary/src/nonlinchem.jl +++ b/lib/ODEProblemLibrary/src/nonlinchem.jl @@ -7,15 +7,15 @@ y0 = [1.0; 0.0; 0.0] tspan = (0.0, 20.0) function nlc_analytic(u0, p, t) [exp(-t); - (2sqrt(exp(-t))besselk(1, 2sqrt(exp(-t))) - - 2besselk(1, 2) / besseli(1, 2) * sqrt(exp(-t))besseli(1, 2sqrt(exp(-t)))) / - (2besselk(0, 2sqrt(exp(-t))) + - (2besselk(1, 2) / besseli(1, 2))besseli(0, 2sqrt(exp(-t)))) - -exp(-t) + 1 + - (-2sqrt(exp(-t)) * besselk(1, 2sqrt(exp(-t))) + - sqrt(exp(-t)) * besseli(1, 2sqrt(exp(-t))) * 2besselk(1, 2) / besseli(1, 2)) / - (2besselk(0, 2sqrt(exp(-t))) + - 2besselk(1, 2) / besseli(1, 2) * besseli(0, 2sqrt(exp(-t))))] + (2sqrt(exp(-t))besselk(1, 2sqrt(exp(-t))) - + 2besselk(1, 2) / besseli(1, 2) * sqrt(exp(-t))besseli(1, 2sqrt(exp(-t)))) / + (2besselk(0, 2sqrt(exp(-t))) + + (2besselk(1, 2) / besseli(1, 2))besseli(0, 2sqrt(exp(-t)))) + -exp(-t) + 1 + + (-2sqrt(exp(-t)) * besselk(1, 2sqrt(exp(-t))) + + sqrt(exp(-t)) * besseli(1, 2sqrt(exp(-t))) * 2besselk(1, 2) / besseli(1, 2)) / + (2besselk(0, 2sqrt(exp(-t))) + + 2besselk(1, 2) / besseli(1, 2) * besseli(0, 2sqrt(exp(-t))))] end nonLinChem_f = ODEFunction(nonLinChem, analytic = nlc_analytic) diff --git a/lib/ODEProblemLibrary/src/ode_linear_prob.jl b/lib/ODEProblemLibrary/src/ode_linear_prob.jl index 4e6055c..d7c8efd 100644 --- a/lib/ODEProblemLibrary/src/ode_linear_prob.jl +++ b/lib/ODEProblemLibrary/src/ode_linear_prob.jl @@ -17,7 +17,7 @@ u(t) = u_0e^{αt} with Float64s. The parameter is ``α`` """ prob_ode_linear = ODEProblem(ODEFunction(linear, analytic = linear_analytic), - 1 / 2, (0.0, 1.0), 1.01) + 1 / 2, (0.0, 1.0), 1.01) """ Linear ODE @@ -35,7 +35,7 @@ u(t) = u_0e^{αt} with BigFloats """ prob_ode_bigfloatlinear = ODEProblem(ODEFunction(linear, analytic = linear_analytic), - big(0.5), (0.0, 1.0), big(1.01)) + big(0.5), (0.0, 1.0), big(1.01)) f_2dlinear = (du, u, p, t) -> (@. du = p * u) f_2dlinear_analytic = (u0, p, t) -> @. u0 * exp(p * t) @@ -56,7 +56,7 @@ u(t) = u_0e^{αt} with Float64s """ prob_ode_2Dlinear = ODEProblem(ODEFunction(f_2dlinear, analytic = f_2dlinear_analytic), - rand(4, 2), (0.0, 1.0), 1.01) + rand(4, 2), (0.0, 1.0), 1.01) """ 100x100 version of the Linear ODE @@ -75,7 +75,7 @@ u(t) = u_0e^{αt} with Float64s """ prob_ode_large2Dlinear = ODEProblem(ODEFunction(f_2dlinear, analytic = f_2dlinear_analytic), - rand(100, 100), (0.0, 1.0), 1.01) + rand(100, 100), (0.0, 1.0), 1.01) """ 4x2 version of the Linear ODE @@ -94,9 +94,9 @@ u(t) = u_0e^{αt} with BigFloats """ prob_ode_bigfloat2Dlinear = ODEProblem(ODEFunction(f_2dlinear, - analytic = f_2dlinear_analytic), - BigFloat.(rand(4, 2)) .* ones(4, 2) / 2, (0.0, 1.0), - big(1.01)) + analytic = f_2dlinear_analytic), + BigFloat.(rand(4, 2)) .* ones(4, 2) / 2, (0.0, 1.0), + big(1.01)) f_2dlinear_notinplace = (u, p, t) -> p * u """ @@ -116,5 +116,5 @@ u(t) = u_0e^{αt} on Float64. Purposefully not in-place as a test. """ prob_ode_2Dlinear_notinplace = ODEProblem(ODEFunction(f_2dlinear_notinplace, - analytic = f_2dlinear_analytic), - rand(4, 2), (0.0, 1.0), 1.01) + analytic = f_2dlinear_analytic), + rand(4, 2), (0.0, 1.0), 1.01) diff --git a/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl b/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl index 6c377fd..a969d83 100644 --- a/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl +++ b/lib/ODEProblemLibrary/src/ode_simple_nonlinear_prob.jl @@ -46,7 +46,7 @@ Fitzhugh-Nagumo (Non-stiff) with initial condition ``v=w=1`` """ prob_ode_fitzhughnagumo = ODEProblem(fitz, [1.0; 1.0], (0.0, 1.0), - (0.7, 0.8, 1 / 12.5, 0.5)) + (0.7, 0.8, 1 / 12.5, 0.5)) #Van der Pol Equations @parameters t μ @@ -166,8 +166,8 @@ Usually solved on ``t₀ = 0.0`` and ``T = 17.0652165601579625588917206249`` Periodic with that setup. """ prob_ode_threebody = ODEProblem(threebody, - [0.994, 0.0, 0.0, big(-2.00158510637908252240537862224)], - (big(0.0), big(17.0652165601579625588917206249))) + [0.994, 0.0, 0.0, big(-2.00158510637908252240537862224)], + (big(0.0), big(17.0652165601579625588917206249))) # Rigid Body Equations @@ -307,36 +307,36 @@ From Hairer Norsett Wanner Solving Ordinary Differential Equations I - Nonstiff Usually solved from 0 to 3. """ prob_ode_pleiades = ODEProblem(pleiades, - [ - 3.0, - 3.0, - -1.0, - -3.0, - 2.0, - -2.0, - 2.0, - 3.0, - -3.0, - 2.0, - 0, - 0, - -4.0, - 4.0, - 0, - 0, - 0, - 0, - 0, - 1.75, - -1.5, - 0, - 0, - 0, - -1.25, - 1, - 0, - 0, - ], (0.0, 3.0)) + [ + 3.0, + 3.0, + -1.0, + -3.0, + 2.0, + -2.0, + 2.0, + 3.0, + -3.0, + 2.0, + 0, + 0, + -4.0, + 4.0, + 0, + 0, + 0, + 0, + 0, + 1.75, + -1.5, + 0, + 0, + 0, + -1.25, + 1, + 0, + 0, + ], (0.0, 3.0)) Random.seed!(100) const mm_A = rand(4, 4) @@ -345,7 +345,7 @@ mm_linear = function (du, u, p, t) end const MM_linear = Matrix(Diagonal(0.5ones(4))) mm_f = ODEFunction(mm_linear; analytic = (u0, p, t) -> exp(inv(MM_linear) * mm_A * t) * u0, - mass_matrix = MM_linear) + mass_matrix = MM_linear) prob_ode_mm_linear = ODEProblem(mm_f, rand(4), (0.0, 1.0)) @parameters t p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 @@ -362,8 +362,8 @@ eqs = [D(y1) ~ -p1 * y1 + p2 * y2 + p3 * y3 + p4, D(y8) ~ -p10 * y6 * y8 + p12 * y7] de = ODESystem(eqs; name = :hires) hires = ODEFunction(de, [y1, y2, y3, y4, y5, y6, y7, y8], - [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], - jac = true) + [p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], + jac = true) u0 = zeros(8) u0[1] = 1 @@ -392,9 +392,9 @@ Reference: [demohires.pdf](http://www.radford.edu/~thompson/vodef90web/problems/ Notebook: [Hires.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Hires.ipynb) """ prob_ode_hires = ODEProblem(hires, u0, (0.0, 321.8122), - (1.71, 0.43, 8.32, 0.0007, 8.75, - 10.03, 0.035, 1.12, 1.745, 280.0, - 0.69, 1.81)) + (1.71, 0.43, 8.32, 0.0007, 8.75, + 10.03, 0.035, 1.12, 1.745, 280.0, + 0.69, 1.81)) @parameters t p1 p2 p3 @variables y1(t) y2(t) y3(t) diff --git a/lib/ODEProblemLibrary/src/pollution_prob.jl b/lib/ODEProblemLibrary/src/pollution_prob.jl index 5473ebd..94d8efe 100644 --- a/lib/ODEProblemLibrary/src/pollution_prob.jl +++ b/lib/ODEProblemLibrary/src/pollution_prob.jl @@ -219,4 +219,4 @@ Reference: [pollu.pdf](https://archimede.dm.uniba.it/~testset/report/pollu.pdf) Notebook: [Pollution.ipynb](http://nbviewer.jupyter.org/github/JuliaDiffEq/DiffEqBenchmarks.jl/blob/master/StiffODE/Pollution.ipynb) """ prob_ode_pollution = ODEProblem(ODEFunction(pollution, jac = pollution_jac), - u0, (0.0, 60.0)) + u0, (0.0, 60.0)) diff --git a/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl b/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl index 87c266a..655a07d 100644 --- a/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl +++ b/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl @@ -9,9 +9,9 @@ RuntimeGeneratedFunctions.init(@__MODULE__) #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, prob_sde_bistable, - prob_sde_bruss, prob_sde_oscilreact + 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 @@ -33,13 +33,13 @@ u(u_0,p,t,W_t)=u_0\exp((α-\frac{β^2}{2})t+βW_t) """ prob_sde_linear = SDEProblem(SDEFunction(f_linear, σ_linear, - analytic = linear_analytic), σ_linear, 1 / 2, - (0.0, 1.0)) + analytic = linear_analytic), σ_linear, 1 / 2, + (0.0, 1.0)) linear_analytic_strat(u0, p, t, W) = @.(u0*exp(1.01t + 0.87W)) prob_sde_linear_stratonovich = SDEProblem(SDEFunction(f_linear, σ_linear, - analytic = linear_analytic_strat), - σ_linear, 1 / 2, (0.0, 1.0)) + analytic = linear_analytic_strat), + σ_linear, 1 / 2, (0.0, 1.0)) f_linear_iip(du, u, p, t) = @.(du=1.01 * u) σ_linear_iip(du, u, p, t) = @.(du=0.87 * u) @doc doc""" @@ -55,11 +55,11 @@ u(u_0,p,t,W_t)=u_0\exp((α-\frac{β^2}{2})t+βW_t) ``` """ prob_sde_2Dlinear = SDEProblem(SDEFunction(f_linear_iip, σ_linear_iip, - analytic = linear_analytic), - σ_linear_iip, ones(4, 2) / 2, (0.0, 1.0)) + analytic = linear_analytic), + σ_linear_iip, ones(4, 2) / 2, (0.0, 1.0)) prob_sde_2Dlinear_stratonovich = SDEProblem(SDEFunction(f_linear_iip, σ_linear_iip, - analytic = linear_analytic_strat), - σ_linear_iip, ones(4, 2) / 2, (0.0, 1.0)) + analytic = linear_analytic_strat), + σ_linear_iip, ones(4, 2) / 2, (0.0, 1.0)) f_cubic(u, p, t) = -0.25 * u * (1 - u^2) σ_cubic(u, p, t) = 0.5 * (1 - u^2) @@ -124,7 +124,7 @@ A multiple dimension extension of `additiveSDEExample` """ prob_sde_additivesystem = SDEProblem(ff_additive_iip, σ_additive_iip, [1.0; 1.0; 1.0; 1.0], - (0.0, 1.0), p) + (0.0, 1.0), p) function f_lorenz(du, u, p, t) du[1] = p[1] * (u[2] - u[1]) @@ -351,9 +351,9 @@ function oval2ModelExample(; largeFluctuations = false, useBigs = false, noiseLe σ = σ2 end - u0 = [0.128483; 1.256853; 0.0030203; 0.0027977; 0.0101511; 0.0422942; 0.2391346; - 0.0008014; 0.0001464; 2.67e-05; 4.8e-6; 9e-7; 0.0619917; 1.2444292; 0.0486676; - 199.9383546; 137.4267984; 1.5180203; 1.5180203] #Fig 9B + u0 = [0.128483 1.256853 0.0030203 0.0027977 0.0101511 0.0422942 0.2391346 + 0.0008014 0.0001464 2.67e-05 4.8e-6 9e-7 0.0619917 1.2444292 0.0486676 + 199.9383546 137.4267984 1.5180203 1.5180203] #Fig 9B if useBigs u0 = big.(u0) end @@ -374,7 +374,7 @@ function stiff_quad_f_ito_analytic(u0, p, t, W) end ff_stiff_quad_ito = SDEFunction(stiff_quad_f_ito, stiff_quad_g, - analytic = stiff_quad_f_ito_analytic) + analytic = stiff_quad_f_ito_analytic) function stiff_quad_f_strat_analytic(u0, p, t, W) α = p[1] @@ -385,7 +385,7 @@ function stiff_quad_f_strat_analytic(u0, p, t, W) end ff_stiff_quad_strat = SDEFunction(stiff_quad_f_strat, stiff_quad_g, - analytic = stiff_quad_f_strat_analytic) + analytic = stiff_quad_f_strat_analytic) @doc doc""" The composite Euler method for stiff stochastic @@ -405,7 +405,7 @@ Higher α or β is stiff, with α being deterministic stiffness and β being noise stiffness (and grows by square). """ prob_sde_stiffquadito = SDEProblem(ff_stiff_quad_ito, stiff_quad_g, 0.5, (0.0, 3.0), - (1.0, 1.0)) + (1.0, 1.0)) @doc doc""" The composite Euler method for stiff stochastic @@ -425,7 +425,7 @@ Higher α or β is stiff, with α being deterministic stiffness and β being noise stiffness (and grows by square). """ prob_sde_stiffquadstrat = SDEProblem(ff_stiff_quad_strat, stiff_quad_g, 0.5, (0.0, 3.0), - (1.0, 1.0)) + (1.0, 1.0)) @doc doc""" Stochastic Heat Equation with scalar multiplicative noise @@ -439,7 +439,7 @@ Raising D or k increases stiffness """ function generate_stiff_stoch_heat(D = 1, k = 1; N = 100) A = full(Tridiagonal([1.0 for i in 1:(N - 1)], [-2.0 for i in 1:N], - [1.0 for i in 1:(N - 1)])) + [1.0 for i in 1:(N - 1)])) dx = 1 / N A = D / (dx^2) * A function f(du, u, p, t) @@ -454,7 +454,7 @@ function generate_stiff_stoch_heat(D = 1, k = 1; N = 100) @. du = k * u end SDEProblem(f, g, ones(N), (0.0, 3.0), - noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3))) + 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] @@ -467,7 +467,7 @@ 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], (0.0, 300.0), p, - noise_rate_prototype = zeros(1, 2)) + 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] @@ -488,7 +488,7 @@ p = (1.0, 1.0, 2.5, 1.0) Stochastic Brusselator """ prob_sde_bruss = SDEProblem(bruss_f, bruss_g, [3.0, 2.0], (0.0, 100.0), p, - noise_rate_prototype = zeros(2, 4)) + noise_rate_prototype = zeros(2, 4)) network = @reaction_network begin p1, (X, Y, Z) --> 0 @@ -508,6 +508,6 @@ p = (0.01, 3.0, 3.0, 4.5, 2.0, 15.0, 20.0, 0.005, 0.01, 0.05) An oscillatory chemical reaction system """ prob_sde_oscilreact = SDEProblem(network, [200.0, 60.0, 120.0, 100.0, 50.0, 50.0, 50.0], - (0.0, 4000.0), p, eval_module = @__MODULE__) + (0.0, 4000.0), p, eval_module = @__MODULE__) end # module