diff --git a/src/ODE.jl b/src/ODE.jl index bd89b27a..d8b80572 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -30,25 +30,25 @@ function hinit(F, x0, t0, p, reltol, abstol) f0 = F(t0, x0) d1 = norm(f0, Inf)/tau - if d0 < 1e-5 || d1 < 1e-5 + if d0 < 1e-5 || d1 < 1e-5 h0 = 1e-6 - else + else h0 = 0.01*d0/d1 - end + end # perform Euler step x1 = x0 + h0*f0 f1 = F(t0 + h0, x1) - # estimate second derivative + # estimate second derivative d2 = norm(f1 - f0, Inf)/(tau*h0) - if max(d1, d2) <= 1e-15 + if max(d1, d2) <= 1e-15 h1 = max(1e-6, 1e-3*h0) - else + else pow = -(2. + log10(max(d1, d2)))/(p + 1.) h1 = 10.^pow - end + end return min(100*h0, h1), f0 end @@ -58,31 +58,31 @@ end # NOTE: in the near future ode23 will be replaced by ode23_bs -#ODE23 Solve non-stiff differential equations. +#ODE23 Solve non-stiff differential equations. # -# ODE23(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system -# of differential equations dy/dt = f(t,y) from t = T0 to t = TFINAL. -# The initial condition is y(T0) = Y0. +# ODE23(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system +# of differential equations dy/dt = f(t,y) from t = T0 to t = TFINAL. +# The initial condition is y(T0) = Y0. # -# The first argument, F, is a function handle or an anonymous function -# that defines f(t,y). This function must have two input arguments, -# t and y, and must return a column vector of the derivatives, dy/dt. +# The first argument, F, is a function handle or an anonymous function +# that defines f(t,y). This function must have two input arguments, +# t and y, and must return a column vector of the derivatives, dy/dt. # -# With two output arguments, [T,Y] = ODE23(...) returns a column -# vector T and an array Y where Y(:,k) is the solution at T(k). +# With two output arguments, [T,Y] = ODE23(...) returns a column +# vector T and an array Y where Y(:,k) is the solution at T(k). # -# More than four input arguments, ODE23(F,TSPAN,Y0,RTOL,P1,P2,...), -# are passed on to F, F(T,Y,P1,P2,...). +# More than four input arguments, ODE23(F,TSPAN,Y0,RTOL,P1,P2,...), +# are passed on to F, F(T,Y,P1,P2,...). # -# ODE23 uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). +# ODE23 uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). # -# Example -# tspan = [0, 2*pi] -# y0 = [1, 0] -# F = (t, y) -> [0 1; -1 0]*y -# ode23(F, tspan, y0) +# Example +# tspan = [0, 2*pi] +# y0 = [1, 0] +# F = (t, y) -> [0 1; -1 0]*y +# ode23(F, tspan, y0) # -# See also ODE23. +# See also ODE23. # Initialize variables. # Adapted from Cleve Moler's textbook @@ -143,7 +143,7 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) y = ynew push!(tout, t) push!(yout, y) - s1 = s4 # Reuse final function value to start new step + s1 = s4 # Reuse final function value to start new step end # Compute a new step size. @@ -175,24 +175,24 @@ end # ode23 # Dormand-Prince 4(5) pair minimizes the local truncation error in the # 5th-order estimate which is what is used to step forward (local extrapolation.) # Generally it produces more accurate results and costs roughly the same -# computationally. The Dormand-Prince pair is the default. +# computationally. The Dormand-Prince pair is the default. # # This is a 4th-order accurate integrator therefore the local error normally -# expected would be O(h^5). However, because this particular implementation +# expected would be O(h^5). However, because this particular implementation # uses the 5th-order estimate for xout (i.e. local extrapolation) moving # forward with the 5th-order estimate should yield errors on the order of O(h^6). # # The order of the RK method is the order of the local *truncation* error, d. # The local truncation error is defined as the principle error term in the -# portion of the Taylor series expansion that gets dropped. This portion of +# portion of the Taylor series expansion that gets dropped. This portion of # the Taylor series exapansion is within the group of terms that gets multipled -# by h in the solution definition of the general RK method. Therefore, the +# by h in the solution definition of the general RK method. Therefore, the # order-p solution created by the RK method will be roughly accurate to -# within O(h^(p+1)). The difference between two different-order solutions is -# the definition of the "local error," l. This makes the local error, l, as +# within O(h^(p+1)). The difference between two different-order solutions is +# the definition of the "local error," l. This makes the local error, l, as # large as the error in the lower order method, which is the truncation # error, d, times h, resulting in O(h^(p+1)). -# Summary: For an RK method of order-p, +# Summary: For an RK method of order-p, # - the local truncation error is O(h^p) # - the local error (used for stepsize adjustment) is O(h^(p+1)) # @@ -207,15 +207,15 @@ end # ode23 # INPUT: # F - User-defined function # Call: xprime = F(t,x) -# t - Time (scalar). -# x - Solution column-vector. +# t - Time (scalar). +# x - Solution column-vector. # xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt. # tspan - [ tstart, tfinal ] # x0 - Initial value COLUMN-vector. # # OUTPUT: -# tout - Returned integration time points (column-vector). -# xout - Returned solution, one solution column-vector per tout-value. +# tout - Returned integration time points (column-vector). +# xout - Returned solution, one solution column-vector per tout-value. # # Original Octave implementation: # Marc Compere @@ -228,10 +228,10 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, maxstep=abs(tspan[end] - tspan[1])/2.5, initstep=0.) # see p.91 in the Ascher & Petzold reference for more infomation. - pow = 1/p # use the higher order to estimate the next step size + pow = 1/p # use the higher order to estimate the next step size - c = sum(a, 2) # consistency condition - k = Array(typeof(x0), length(c)) # storage for intermediate steps + c = sum(a, 2) # consistency condition + k = Array(typeof(x0), length(c)) # storage for intermediate steps # Initialization t = tspan[1] @@ -240,10 +240,10 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, h = initstep if h == 0. - # initial guess at a step size - h, k[1] = hinit(F, x0, t, p, reltol, abstol) - else - k[1] = F(t,x0) # first stage + # initial guess at a step size + h, k[1] = hinit(F, x0, t, p, reltol, abstol) + else + k[1] = F(t,x0) # first stage end h = tdir*min(h, maxstep) @@ -278,8 +278,8 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8, gamma1 = xs - xp # Estimate the error and the acceptable error - delta = norm(gamma1, Inf) # actual error - tau = max(reltol*norm(x,Inf),abstol) # allowable error + delta = norm(gamma1, Inf) # actual error + tau = max(reltol*norm(x,Inf),abstol) # allowable error # Update the solution only if the error is acceptable if delta <= tau @@ -316,10 +316,10 @@ end # Bogacki–Shampine coefficients const bs_coefficients = (3, - [ 0 0 0 0 - 1/2 0 0 0 - 0 3/4 0 0 - 2/9 1/3 4/9 0], + [ 0 0 0 0 + 1/2 0 0 0 + 0 3/4 0 0 + 2/9 1/3 4/9 0], # 2nd order b-coefficients [7/24 1/4 1/3 1/8], # 3rd order b-coefficients @@ -329,19 +329,19 @@ ode23_bs(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, bs_coefficients...; kwa # Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableau in -# U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations +# U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations # and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics # (SIAM), Philadelphia, 1998 # # Dormand-Prince coefficients const dp_coefficients = (5, - [ 0 0 0 0 0 0 - 1/5 0 0 0 0 0 - 3/40 9/40 0 0 0 0 - 44/45 -56/15 32/9 0 0 0 - 19372/6561 -25360/2187 64448/6561 -212/729 0 0 - 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0 - 35/384 0 500/1113 125/192 -2187/6784 11/84], + [ 0 0 0 0 0 0 + 1/5 0 0 0 0 0 + 3/40 9/40 0 0 0 0 + 44/45 -56/15 32/9 0 0 0 + 19372/6561 -25360/2187 64448/6561 -212/729 0 0 + 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0 + 35/384 0 500/1113 125/192 -2187/6784 11/84], # 4th order b-coefficients [5179/57600 0 7571/16695 393/640 -92097/339200 187/2100 1/40], # 5th order b-coefficients @@ -352,12 +352,12 @@ ode45_dp(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, dp_coefficients...; kwa # Fehlberg coefficients const fb_coefficients = (5, - [ 0 0 0 0 0 - 1/4 0 0 0 0 - 3/32 9/32 0 0 0 - 1932/2197 -7200/2197 7296/2197 0 0 - 439/216 -8 3680/513 -845/4104 0 - -8/27 2 -3544/2565 1859/4104 -11/40], + [ 0 0 0 0 0 + 1/4 0 0 0 0 + 3/32 9/32 0 0 0 + 1932/2197 -7200/2197 7296/2197 0 0 + 439/216 -8 3680/513 -845/4104 0 + -8/27 2 -3544/2565 1859/4104 -11/40], # 4th order b-coefficients [25/216 0 1408/2565 2197/4104 -1/5 0], # 5th order b-coefficients @@ -369,12 +369,12 @@ ode45_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients...; kwa # Cash-Karp coefficients # Numerical Recipes in Fortran 77 const ck_coefficients = (5, - [ 0 0 0 0 0 - 1/5 0 0 0 0 - 3/40 9/40 0 0 0 - 3/10 -9/10 6/5 0 0 - -11/54 5/2 -70/27 35/27 0 - 1631/55296 175/512 575/13824 44275/110592 253/4096], + [ 0 0 0 0 0 + 1/5 0 0 0 0 + 3/40 9/40 0 0 0 + 3/10 -9/10 6/5 0 0 + -11/54 5/2 -70/27 35/27 0 + 1631/55296 175/512 575/13824 44275/110592 253/4096], # 4th order b-coefficients [37/378 0 250/621 125/594 0 512/1771], # 5th order b-coefficients @@ -387,19 +387,19 @@ ode45_ck(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, ck_coefficients...; kwa # Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control". # National Aeronautics and Space Administration. const fb_coefficients_78 = (8, - [ 0 0 0 0 0 0 0 0 0 0 0 0 - 2/27 0 0 0 0 0 0 0 0 0 0 0 - 1/36 1/12 0 0 0 0 0 0 0 0 0 0 - 1/24 0 1/8 0 0 0 0 0 0 0 0 0 - 5/12 0 -25/16 25/16 0 0 0 0 0 0 0 0 - 1/20 0 0 1/4 1/5 0 0 0 0 0 0 0 - -25/108 0 0 125/108 -65/27 125/54 0 0 0 0 0 0 - 31/300 0 0 0 61/225 -2/9 13/900 0 0 0 0 0 - 2 0 0 -53/6 704/45 -107/9 67/90 3 0 0 0 0 - -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0 0 0 - 2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 - 3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 - -1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 1], + [ 0 0 0 0 0 0 0 0 0 0 0 0 + 2/27 0 0 0 0 0 0 0 0 0 0 0 + 1/36 1/12 0 0 0 0 0 0 0 0 0 0 + 1/24 0 1/8 0 0 0 0 0 0 0 0 0 + 5/12 0 -25/16 25/16 0 0 0 0 0 0 0 0 + 1/20 0 0 1/4 1/5 0 0 0 0 0 0 0 + -25/108 0 0 125/108 -65/27 125/54 0 0 0 0 0 0 + 31/300 0 0 0 61/225 -2/9 13/900 0 0 0 0 0 + 2 0 0 -53/6 704/45 -107/9 67/90 3 0 0 0 0 + -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0 0 0 + 2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 + 3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 + -1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 1], # 7th order b-coefficients [41/840 0 0 0 0 34/105 9/35 9/35 9/280 9/280 41/840 0 0], # 8th order b-coefficients @@ -417,15 +417,15 @@ const ode45 = ode45_dp # P.J. Prince and J.R.Dormand: High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 1981. -#ODE4 Solve non-stiff differential equations, fourth order -# fixed-step Runge-Kutta method. +#ODE4 Solve non-stiff differential equations, fourth order +# fixed-step Runge-Kutta method. # -# [T,X] = ODE4(ODEFUN, X0, TSPAN) with TSPAN = [T0:H:TFINAL] -# integrates the system of differential equations x' = f(t,x) from time -# T0 to TFINAL in steps of H with initial conditions X0. Function -# ODEFUN(T,X) must return a column vector corresponding to f(t,x). Each -# row in the solution array X corresponds to a time returned in the -# column vector T. +# [T,X] = ODE4(ODEFUN, X0, TSPAN) with TSPAN = [T0:H:TFINAL] +# integrates the system of differential equations x' = f(t,x) from time +# T0 to TFINAL in steps of H with initial conditions X0. Function +# ODEFUN(T,X) must return a column vector corresponding to f(t,x). Each +# row in the solution array X corresponds to a time returned in the +# column vector T. function ode4(F, x0, tspan) h = diff(tspan) x = Array(typeof(x0), length(tspan)) @@ -475,7 +475,7 @@ function fdjacobian(F, x, t) return dFdx end -# ODE23S Solve stiff systems based on a modified Rosenbrock triple +# ODE23S Solve stiff systems based on a modified Rosenbrock triple # (also used by MATLAB's ODE23s); see Sec. 4.1 in # # [SR97] L.F. Shampine and M.W. Reichelt: "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22 @@ -511,10 +511,10 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, h = initstep if h == 0. - # initial guess at a step size - h, F0 = hinit(F, y0, t, 3, reltol, abstol) - else - F0 = F(t,y) + # initial guess at a step size + h, F0 = hinit(F, y0, t, 3, reltol, abstol) + else + F0 = F(t,y) end h = tdir*min(h, maxstep) @@ -535,7 +535,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, W = one(J) - h*d*J else # note: if there is a mass matrix M on the lhs of the ODE, i.e., - # M * dy/dt = F(t,y) + # M * dy/dt = F(t,y) # we can simply replace eye(J) by M in the following expression # (see Sec. 5 in [SR97]) @@ -557,7 +557,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8, delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error # check if new solution is acceptable - if err <= delta + if err <= delta if points==:specified || points==:all # only points in tspan are requested @@ -595,7 +595,7 @@ end #ODEROSENBROCK Solve stiff differential equations, Rosenbrock method -# with provided coefficients. +# with provided coefficients. function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing) if typeof(jacobian) == Function @@ -643,42 +643,42 @@ end # Kaps-Rentrop coefficients const kr4_coefficients = (0.231, - [0 0 0 0 - 2 0 0 0 - 4.452470820736 4.16352878860 0 0 - 4.452470820736 4.16352878860 0 0], - [3.95750374663 4.62489238836 0.617477263873 1.28261294568], - [ 0 0 0 0 - -5.07167533877 0 0 0 - 6.02015272865 0.1597500684673 0 0 - -1.856343618677 -8.50538085819 -2.08407513602 0],) + [0 0 0 0 + 2 0 0 0 + 4.452470820736 4.16352878860 0 0 + 4.452470820736 4.16352878860 0 0], + [3.95750374663 4.62489238836 0.617477263873 1.28261294568], + [ 0 0 0 0 + -5.07167533877 0 0 0 + 6.02015272865 0.1597500684673 0 0 + -1.856343618677 -8.50538085819 -2.08407513602 0],) ode4s_kr(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, kr4_coefficients...; jacobian=jacobian) # Shampine coefficients const s4_coefficients = (0.5, [ 0 0 0 0 - 2 0 0 0 - 48/25 6/25 0 0 - 48/25 6/25 0 0], + 2 0 0 0 + 48/25 6/25 0 0 + 48/25 6/25 0 0], [19/9 1/2 25/108 125/108], - [ 0 0 0 0 - -8 0 0 0 - 372/25 12/5 0 0 - -112/125 -54/125 -2/5 0],) + [ 0 0 0 0 + -8 0 0 0 + 372/25 12/5 0 0 + -112/125 -54/125 -2/5 0],) ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coefficients...; jacobian=jacobian) # Use Shampine coefficients by default (matching Numerical Recipes) const ode4s = ode4s_s -const ms_coefficients4 = [ 1 0 0 0 - -1/2 3/2 0 0 - 5/12 -4/3 23/12 0 - -9/24 37/24 -59/24 55/24] +const ms_coefficients4 = [ 1 0 0 0 + -1/2 3/2 0 0 + 5/12 -4/3 23/12 0 + -9/24 37/24 -59/24 55/24] # ODE_MS Fixed-step, fixed-order multi-step numerical method -# with Adams-Bashforth-Moulton coefficients +# with Adams-Bashforth-Moulton coefficients function ode_ms(F, x0, tspan, order::Integer) h = diff(tspan) x = Array(typeof(x0), length(tspan)) @@ -692,10 +692,10 @@ function ode_ms(F, x0, tspan, order::Integer) for s = 5:order for j = 0:(s - 1) # Assign in correct order for multiplication below - # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) + # (a factor depending on j and s) .* (an integral of a polynomial with -(0:s), except -j, as roots) p_int = polyint(poly(diagm(-[0:j - 1; j + 1:s - 1]))) b[s, s - j] = ((-1)^j / factorial(j) - / factorial(s - 1 - j) * polyval(p_int, 1)) + / factorial(s - 1 - j) * polyval(p_int, 1)) end end end