Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

AB adaptive #2

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

AB adaptive #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
30dce73
Update Simple_Differential_Equation.ipynb
jpfairbanks Mar 12, 2016
8c2ec0d
On tableau was in terms of 1/2 instead of 1//2, etc.
mauro3 May 30, 2016
e5c46a2
example of a user defined type
Jul 11, 2016
63ad4b0
Updating example with comments and Compat
Jul 13, 2016
75bfc49
Update README.md
ChrisRackauckas Sep 1, 2016
2394aee
Update README.md
ChrisRackauckas Sep 1, 2016
4bd341b
Merge pull request #97 from mauro3/m3/bt_rk23-tableau-fix
ChrisRackauckas Sep 1, 2016
f6fc008
Merge pull request #107 from jkozdon/master
ChrisRackauckas Sep 1, 2016
0d1a3d1
Merge pull request #92 from jpfairbanks/patch-1
ChrisRackauckas Sep 2, 2016
145e081
Addition of Variable and Fixed Adam Moulton Method
obiajulu May 23, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
Various basic Ordinary Differential Equation solvers implemented in Julia.

[![Build Status](https://travis-ci.org/JuliaLang/ODE.jl.svg?branch=master)](https://travis-ci.org/JuliaLang/ODE.jl)
[![Coverage Status](https://img.shields.io/coveralls/JuliaLang/ODE.jl.svg)](https://coveralls.io/r/JuliaLang/ODE.jl)
[![ODE](http://pkg.julialang.org/badges/ODE_0.3.svg)](http://pkg.julialang.org/?pkg=ODE&ver=0.3)
[![ODE](http://pkg.julialang.org/badges/ODE_0.4.svg)](http://pkg.julialang.org/?pkg=ODE&ver=0.4)
[![Join the chat at https://gitter.im/JuliaDiffEq/Lobby](https://badges.gitter.im/JuliaDiffEq/Lobby.svg)](https://gitter.im/JuliaDiffEq/Lobby?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge)
[![Build Status](https://travis-ci.org/JuliaDiffEq/ODE.jl.svg?branch=master)](https://travis-ci.org/JuliaDiffEq/ODE.jl)
[![Coverage Status](https://img.shields.io/coveralls/JuliaDiffEq/ODE.jl.svg)](https://coveralls.io/r/JuliaDiffEq/ODE.jl)
[![ODE](http://pkg.julialang.org/badges/ODE_0.4.svg)](http://pkg.julialang.org/?pkg=ODE)
[![ODE](http://pkg.julialang.org/badges/ODE_0.5.svg)](http://pkg.julialang.org/?pkg=ODE)

Pull requests are always highly welcome to fix bugs, add solvers, or anything else!

Expand Down
4 changes: 2 additions & 2 deletions examples/Simple_Differential_Equation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@
"outputs": [],
"source": [
"# Initial condtions -- x(0) and x'(0)\n",
"const start = [0.0; 0.1]\n",
"const starty = [0.0; 0.1]\n",
"\n",
"# Time vector going from 0 to 2*PI in 0.01 increments\n",
"time = 0:0.1:4*pi;"
Expand All @@ -113,7 +113,7 @@
},
"outputs": [],
"source": [
"t, y = ode45(f, start, time);"
"t, y = ode45(f, starty, time);"
]
},
{
Expand Down
53 changes: 53 additions & 0 deletions examples/User_Type.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
#=
This is an example of how to use a custom user type of with ODE module.

An additional example can be found in `test/inverface-tests.jl`

The ODE considered in this example is the simple oscilator

a' = -b
b' = a

with initial condtion a(t=0) = 0 and b(t=0) = 0.1.
=#

using ODE
using PyPlot
using Compat

typealias DataFloat Float64
type Vec2
a::DataFloat
b::DataFloat
end

# This is the minimal set of function required on the type inorder to work with
# the ODE module
@compat Base.:/(x::Vec2, y::Real) = Vec2(x.a/y, x.b/y)
@compat Base.:*(y::Real, x::Vec2) = Vec2(y * x.a, y * x.b)
@compat Base.:*(x::Vec2, y::Real) = y*x
@compat Base.:.*(y::Real, x::Vec2) = y*x
@compat Base.:+(x::Vec2, y::Vec2) = Vec2(x.a + y.a, x.b + y.b)
@compat Base.:-(x::Vec2, y::Vec2) = Vec2(x.a - y.a, x.b - y.b)
Base.norm(x::Vec2) = sqrt(x.a^2 + x.b^2)
Base.zero(x::Type{Vec2}) = Vec2(zero(DataFloat), zero(DataFloat))
ODE.isoutofdomain(x::Vec2) = isnan(x.a) || isnan(x.b)

# RHS function
f(t,y) = Vec2(-y.b, y.a)

# Initial condtions
start = Vec2(0.0, 0.1)

# Time vector going from 0 to 2*PI in 0.01 increments
time = 0:0.1:4*pi

# Solve the ODE
t, y = ode45(f, start, time)

# Plot the solution
a = map(y -> y.a, y)
b = map(y -> y.b, y)
plot(t, a, label="a(t)")
plot(t, b, label="b(t)")
legend()
68 changes: 17 additions & 51 deletions src/ODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,9 @@ using Compat

## minimal function export list
# adaptive non-stiff:
export ode23, ode45, ode78
export ode23, ode45, ode78, ode113
# non-adaptive non-stiff:
export ode4, ode4ms
export ode4, ode4ms, ode4am
# adaptive stiff:
export ode23s
# non-adaptive stiff:
Expand Down Expand Up @@ -93,7 +93,7 @@ end
# of the allowed domain. Used in adaptive step-control.
isoutofdomain(x) = isnan(x)

function make_consistent_types(fn, y0, tspan, btab::Tableau)
function make_consistent_types(fn, y0, tspan)
# There are a few types involved in a call to a ODE solver which
# somehow need to be consistent:
#
Expand All @@ -113,7 +113,6 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau)
# - Eyf: suitable eltype of y and f(t,y)
# --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))
# - Ty: container type of y0
# - btab: tableau with entries converted to Et

# Needed interface:
# On components: /, -
Expand All @@ -138,6 +137,19 @@ function make_consistent_types(fn, y0, tspan, btab::Tableau)
!isleaftype(Et) && warn("The eltype(tspan) is not a concrete type! Change type of tspan for better performance.")
!isleaftype(Eyf) && warn("The eltype(y0/tspan[1]) is not a concrete type! Change type of y0 and/or tspan for better performance.")

return Et, Eyf, Ty
end


function make_consistent_types(fn, y0, tspan, btab::Tableau)
# Returns
# - Et: eltype of time, needs to be a real "continuous" type, at
# the moment a AbstractFloat
# - Eyf: suitable eltype of y and f(t,y)
# --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))
# - Ty: container type of y0
# - btab: tableau with entries converted to Et
Et, Eyf, Ty = make_consistent_types(fn,y0,tspan)
btab_ = convert(Et, btab)
return Et, Eyf, Ty, btab_
end
Expand All @@ -147,48 +159,7 @@ end
###############################################################################

include("runge_kutta.jl")

# ODE_MS Fixed-step, fixed-order multi-step numerical method
# with Adams-Bashforth-Moulton coefficients
function ode_ms(F, x0, tspan, order::Integer)
h = diff(tspan)
x = Array(typeof(x0), length(tspan))
x[1] = x0

if 1 <= order <= 4
b = ms_coefficients4
else
b = zeros(order, order)
b[1:4, 1:4] = ms_coefficients4
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)
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))
end
end
end

# TODO: use a better data structure here (should be an order-element circ buffer)
xdot = similar(x)
for i = 1:length(tspan)-1
# Need to run the first several steps at reduced order
steporder = min(i, order)
xdot[i] = F(tspan[i], x[i])

x[i+1] = x[i]
for j = 1:steporder
x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)]
end
end
return vcat(tspan), x
end

# Use order 4 by default
ode4ms(F, x0, tspan) = ode_ms(F, x0, tspan, 4)
ode5ms(F, x0, tspan) = ODE.ode_ms(F, x0, tspan, 5)
include("adams_methods.jl")

###############################################################################
## STIFF SOLVERS
Expand Down Expand Up @@ -419,10 +390,5 @@ ode4s_s(F, x0, tspan; jacobian=nothing) = oderosenbrock(F, x0, tspan, s4_coeffic
# 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]


end # module ODE
Loading