From 1dcd171d074757463ef33931610927dcb8fdefac Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Fri, 20 Mar 2015 09:26:54 -0600 Subject: [PATCH 1/6] Added parameter splatting to ode23 --- src/ODE.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 8811574fe..548de92b5 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -48,7 +48,7 @@ export ode4s, ode4ms, ode4 # Initialize variables. # Adapted from Cleve Moler's textbook # http://www.mathworks.com/moler/ncm/ode23tx.m -function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) +function ode23(F, y0, tspan, params...; reltol = 1.e-5, abstol = 1.e-8) if reltol == 0 warn("setting reltol = 0 gives a step size of zero") end @@ -69,7 +69,7 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) # Compute initial step size. - s1 = F(t, y) + s1 = F(t, y, params...) r = norm(s1./max(abs(y), threshold), Inf) + realmin() # TODO: fix type bug in max() h = tdir*0.8*reltol^(1/3)/r @@ -89,11 +89,11 @@ function ode23(F, y0, tspan; reltol = 1.e-5, abstol = 1.e-8) # Attempt a step. - s2 = F(t+h/2, y+h/2*s1) - s3 = F(t+3*h/4, y+3*h/4*s2) + s2 = F(t+h/2, y+h/2*s1, params...) + s3 = F(t+3*h/4, y+3*h/4*s2, params...) tnew = t + h ynew = y + h*(2*s1 + 3*s2 + 4*s3)/9 - s4 = F(tnew, ynew) + s4 = F(tnew, ynew, params...) # Estimate the error. From 26982074b06f7e56ce4e43e9a7c171522db2c426 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Fri, 20 Mar 2015 09:34:47 -0600 Subject: [PATCH 2/6] Updated documentation of ode23 to have correct order of arguments in the example, reduce the MATLAB-ness, add the new parameter functionality, comment that there are keyword arguments --- src/ODE.jl | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 548de92b5..86f4969e9 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -19,38 +19,42 @@ export ode4s, ode4ms, ode4 # ode4ms, ode_ms -#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, y0, tspan) with tspan = [t0, t_final] integrates the system +# of differential equations dy/dt = f(t,y) from t = t0 to t = t_final. +# Here, y may be a scalar or a vector. +# 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. +# t and y, and must return a 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). +# ode23 returns the pair (tout, yout), +# where tout is the vector of times and yout an array of solutions: +# yout[k,:] is the solution at time tout(k) # -# More than four input arguments, ODE23(F,TSPAN,Y0,RTOL,P1,P2,...), -# are passed on to F, F(T,Y,P1,P2,...). +# Parameters may be passed through to the function F by adding additional +# arguments: +# ode23(f, y0, tspan, p1, p2, ...) +# calls f(t, y, p1, p2, ...). # -# ODE23 uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). +# Keyword arguments reltol and abstol specify the relative and absolute +# tolerances, respectively; their default values are reltol=1.e-5; abstol=1.e-8. # -# Example +# 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) -# -# See also ODE23. +# ode23(F, y0, tspan) -# Initialize variables. # Adapted from Cleve Moler's textbook # http://www.mathworks.com/moler/ncm/ode23tx.m function ode23(F, y0, tspan, params...; reltol = 1.e-5, abstol = 1.e-8) if reltol == 0 - warn("setting reltol = 0 gives a step size of zero") + warn("Setting reltol = 0 gives a step size of zero") end threshold = abstol / reltol From e7591c98d0842e0103b419e75b582e34d5a9e053 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Mon, 23 Mar 2015 10:08:42 -0600 Subject: [PATCH 3/6] Added test of ode23 with parameters --- test/runtests.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 7e86984c5..b970d1479 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,6 +52,19 @@ for solver in solvers @test maximum(abs(ys-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol end +# Test ode23 with parameters: + +solver = ode23 +println("using $solver with parameters") +# dy +# -- = -αy ==> y = y0*e.^(-αt) +# dt + +for α in 1.0:1.0:10. + t,y=solver((t,y,α)->-α*y, 1., [0:.001:1;], α) + @test maximum(abs(y-e.^(-α*t))) < tol +end + # rober testcase from http://www.unige.ch/~hairer/testset/testset.html let println("ROBER test case") @@ -71,4 +84,7 @@ let @test norm(refsol-y[end], Inf) < 2e-10 end + + + println("All looks OK") From 89ea090a314b649f655240332bc4f415894ff875 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Mon, 23 Mar 2015 10:16:28 -0600 Subject: [PATCH 4/6] Added example of passing parameters to documentation. Added '.' to numerical values so that the example code does not give InexactErrors. Removed spaces in docs of ode23 --- src/ODE.jl | 54 +++++++++++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 86f4969e9..05765c48a 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -13,7 +13,7 @@ export ode23s export ode4s, ode4ms, ode4 ## complete function export list -#export ode23, ode4, +# export ode23, ode4, # oderkf, ode45, ode45_dp, ode45_fb, ode45_ck, # oderosenbrock, ode4s, ode4s_kr, ode4s_s, # ode4ms, ode_ms @@ -21,37 +21,45 @@ export ode4s, ode4ms, ode4 # ode23: Solve non-stiff differential equations. # -# ode23(f, y0, tspan) with tspan = [t0, t_final] integrates the system -# of differential equations dy/dt = f(t,y) from t = t0 to t = t_final. -# Here, y may be a scalar or a vector. -# The initial condition is y(t0) = y0. +# ode23(f, y0, tspan) with tspan = [t0, t_final] integrates the system +# of differential equations dy/dt = f(t,y) from t = t0 to t = t_final. +# Here, y may be a scalar or a vector. +# 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 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 vector of the derivatives, dy/dt. # -# ode23 returns the pair (tout, yout), -# where tout is the vector of times and yout an array of solutions: -# yout[k,:] is the solution at time tout(k) +# ode23 returns the pair (tout, yout), +# where tout is the vector of times and yout an array of solutions: +# yout[k,:] is the solution at time tout(k) # -# Parameters may be passed through to the function F by adding additional -# arguments: -# ode23(f, y0, tspan, p1, p2, ...) -# calls f(t, y, p1, p2, ...). +# Parameters may be passed through to the function F by adding additional +# arguments: +# ode23(f, y0, tspan, p1, p2, ...) +# calls f(t, y, p1, p2, ...). # -# Keyword arguments reltol and abstol specify the relative and absolute -# tolerances, respectively; their default values are reltol=1.e-5; abstol=1.e-8. +# Keyword arguments reltol and abstol specify the relative and absolute +# tolerances, respectively; their default values are reltol=1.e-5; abstol=1.e-8. # -# 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, y0, tspan) +# Example: +# tspan = [0., 2*pi] +# y0 = [1.0, 0.0] +# F = (t, y) -> [0.0 1.0; -1.0 0.0]*y +# ode23(F, y0, tspan) +# +# Example of passing through parameters: +# tspan = [0, 2*pi] +# F(t, y, α, β) = -α*y + β +# y0 = 1.0 +# α, β = 1.0, 0.5 +# ode23(F, y0, tspan, α, β) # Adapted from Cleve Moler's textbook # http://www.mathworks.com/moler/ncm/ode23tx.m + function ode23(F, y0, tspan, params...; reltol = 1.e-5, abstol = 1.e-8) if reltol == 0 warn("Setting reltol = 0 gives a step size of zero") From 0980c623a3324dd4f480e6de88c84a4e671a47f8 Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Mon, 23 Mar 2015 10:58:26 -0600 Subject: [PATCH 5/6] Converted the documentation of ode23 to use @doc --- src/ODE.jl | 90 +++++++++++++++++++++++++++++------------------------- 1 file changed, 49 insertions(+), 41 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 05765c48a..8e9a713b3 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -2,6 +2,8 @@ module ODE +VERSION < v"0.4-" && using Docile + using Polynomials ## minimal function export list @@ -19,47 +21,53 @@ export ode4s, ode4ms, ode4 # ode4ms, ode_ms -# ode23: Solve non-stiff differential equations. -# -# ode23(f, y0, tspan) with tspan = [t0, t_final] integrates the system -# of differential equations dy/dt = f(t,y) from t = t0 to t = t_final. -# Here, y may be a scalar or a vector. -# 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 vector of the derivatives, dy/dt. -# -# ode23 returns the pair (tout, yout), -# where tout is the vector of times and yout an array of solutions: -# yout[k,:] is the solution at time tout(k) -# -# Parameters may be passed through to the function F by adding additional -# arguments: -# ode23(f, y0, tspan, p1, p2, ...) -# calls f(t, y, p1, p2, ...). -# -# Keyword arguments reltol and abstol specify the relative and absolute -# tolerances, respectively; their default values are reltol=1.e-5; abstol=1.e-8. -# -# ode23 uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). -# -# Example: -# tspan = [0., 2*pi] -# y0 = [1.0, 0.0] -# F = (t, y) -> [0.0 1.0; -1.0 0.0]*y -# ode23(F, y0, tspan) -# -# Example of passing through parameters: -# tspan = [0, 2*pi] -# F(t, y, α, β) = -α*y + β -# y0 = 1.0 -# α, β = 1.0, 0.5 -# ode23(F, y0, tspan, α, β) - -# Adapted from Cleve Moler's textbook -# http://www.mathworks.com/moler/ncm/ode23tx.m - +@doc doc""" +`ode23`: Solve non-stiff differential equations. + +`ode23(f, y0, tspan)` with `tspan = [t0, t_final]` integrates the system +of differential equations \$dy/dt = f(t,y)\$ from \$t = t_0\$ to \$t = t_\mathrm{final}\$. +Here, \$y\$ may be a scalar or a vector. +The initial condition is \$y(t0) = y_0\$. + +The first argument, `F`, is a function that defines \$f(t,y)\$. +This function must have at least two input arguments, +`t` and `y`, and must return a vector of the derivatives, \$dy_i/dt\$. + +`ode23` returns the tuple `(tout, yout)`, +where `tout` is the vector of times and yout an array of solutions: +`yout[k,:]` is the solution at time `tout(k)`. + +Parameters may be passed through to the function `F` by adding additional +arguments: +`ode23(f, y0, tspan, p1, p2, ...)` +calls `f(t, y, p1, p2, ...)`. + +Keyword arguments `reltol` and `abstol` specify the relative and absolute +tolerances, respectively; their default values are `reltol=1.e-5` +and `abstol=1.e-8`. + +`ode23` uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23). + +Example of the basic use of `ode23`: +``` + tspan = [0., 2*pi] + y0 = [1.0, 0.0] + F = (t, y) -> [0.0 1.0; -1.0 0.0]*y + ode23(F, y0, tspan) +``` + +Example of passing through parameters: +``` + tspan = [0, 2*pi] + F(t, y, α, β) = -α*y + β + y0 = 1.0 + α, β = 1.0, 0.5 + ode23(F, y0, tspan, α, β) +``` + +The code for `ode23` was adapted from Cleve Moler's textbook: + +""" -> function ode23(F, y0, tspan, params...; reltol = 1.e-5, abstol = 1.e-8) if reltol == 0 warn("Setting reltol = 0 gives a step size of zero") From d991e1d4922742be83e0afc5392eccf60eb2d61d Mon Sep 17 00:00:00 2001 From: "David P. Sanders" Date: Mon, 23 Mar 2015 11:30:49 -0600 Subject: [PATCH 6/6] Added Docile to REQUIRE --- REQUIRE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/REQUIRE b/REQUIRE index a9aa1a122..cd1a90fad 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,4 @@ julia 0.2 Polynomials +Docile +