Skip to content
This repository has been archived by the owner on Sep 9, 2024. It is now read-only.

Commit

Permalink
fix for v1.0
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Oct 2, 2018
1 parent 0f1aaf6 commit 68227a1
Show file tree
Hide file tree
Showing 7 changed files with 22 additions and 28 deletions.
2 changes: 1 addition & 1 deletion REQUIRE
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
julia 0.7-beta2
julia 0.7
Polynomials
Compat 0.17.0
DiffEqBase 4.0.0
Expand Down
6 changes: 3 additions & 3 deletions src/ODE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ function ode_ms(F, x0, tspan, order::Integer; kwargs...)
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])))
p_int = polyint(poly(diagm(0 => -[0:j - 1; j + 1:s - 1])))
b[s, s - j] = ((-1)^j / factorial(j)
/ factorial(s - 1 - j) * polyval(p_int, 1))
end
Expand Down Expand Up @@ -254,7 +254,7 @@ function ode23s(F, y0, tspan;
reltol = 1.0e-5, abstol = 1.0e-8,
jacobian=nothing,
points=:all,
norm=Base.norm,
norm=LinearAlgebra.norm,
minstep=abs(tspan[end] - tspan[1])/1e18,
maxstep=abs(tspan[end] - tspan[1])/2.5,
initstep=0.)
Expand Down Expand Up @@ -307,7 +307,7 @@ function ode23s(F, y0, tspan;
# we can simply replace eye(J) by M in the following expression
# (see Sec. 5 in [SR97])

W = lufact( I - h*d*J )
W = lu( I - h*d*J )
end

# approximate time-derivative of F
Expand Down
2 changes: 1 addition & 1 deletion src/common.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ function solve(
p = prob.p

if isinplace
f = (t,u) -> (du = zeros(u); prob.f(du,u,p,t); vec(du))
f = (t,u) -> (du = zero(u); prob.f(du,u,p,t); vec(du))
elseif uType <: AbstractArray
f = (t,u) -> vec(prob.f(reshape(u,sizeu),p,t))
else
Expand Down
4 changes: 2 additions & 2 deletions src/runge_kutta.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ conv_field(D,a::Array{T,N}) where {T,N} = convert(Array{D,N}, a)
function Base.convert(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T}) where {Tnew<:Real,Name,S,T}
# Converts the tableau coefficients to the new type Tnew
newflds = ()
for n in fieldnames(tab)
for n in fieldnames(typeof(tab))
fld = getfield(tab,n)
if eltype(fld)==T
newflds = tuple(newflds..., conv_field(Tnew, fld))
Expand Down Expand Up @@ -227,7 +227,7 @@ function oderk_adapt(fn, y0, tspan, btab::TableauRKExplicit; kwords...)
end
function oderk_adapt(fn, y0::AbstractVector, tspan, btab_::TableauRKExplicit{N,S};
reltol = 1.0e-5, abstol = 1.0e-8,
norm=Base.norm,
norm=LinearAlgebra.norm,
minstep=abs(tspan[end] - tspan[1])/1e18,
maxstep=abs(tspan[end] - tspan[1])/2.5,
initstep=0,
Expand Down
2 changes: 1 addition & 1 deletion test/common.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using ODE, DiffEqBase, DiffEqProblemLibrary
using ODE, DiffEqBase, DiffEqProblemLibrary, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
using LinearAlgebra
Expand Down
24 changes: 9 additions & 15 deletions test/interface-tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,26 +4,29 @@
# This is to test a scalar-like state variable
# (due to @acroy: https://gist.github.com/acroy/28be4f2384d01f38e577)

import Base: +, -, *, /, .+, .-, .*, ./
import Base: +, -, *, /
using LinearAlgebra

const delta0 = 0.
const V0 = 1.
const g0 = 0.

# define custom type ...
struct CompSol
struct CompSol <: AbstractMatrix{ComplexF64}
rho::Matrix{ComplexF64}
x::Float64
p::Float64

CompSol(r,x,p) = new(copy(r),x,p)
end

Base.getindex(A::CompSol,i...) = A.rho[i...]
Base.setindex!(A,x,i...) = setindex!(A.rho,x,i...)

# ... which has to support the following operations
# to work with odeX
Base.norm(y::CompSol, p::Float64) = maximum([Base.norm(y.rho, p) abs(y.x) abs(y.p)])
Base.norm(y::CompSol) = norm(y::CompSol, 2.0)
LinearAlgebra.norm(y::CompSol, p::Float64) = maximum([LinearAlgebra.norm(y.rho, p) abs(y.x) abs(y.p)])
LinearAlgebra.norm(y::CompSol) = LinearAlgebra.norm(y::CompSol, 2.0)

+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
Expand All @@ -32,19 +35,10 @@ Base.norm(y::CompSol) = norm(y::CompSol, 2.0)
/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)

### new for PR #68
Base.abs(y::CompSol) = norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia
Base.abs(y::CompSol) = LinearAlgebra.norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia
Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.)
ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p))

# Because the new RK solvers wrap scalars in an array and because of
# https://github.com/JuliaLang/julia/issues/11053 these are also needed:
.+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
.-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
.*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s)
.*(s::Real, y1::CompSol) = y1*s
./(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)


################################################################################

# define RHSs of differential equations
Expand Down Expand Up @@ -75,7 +69,7 @@ for solver in solvers
continue
end
t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, range(0., stop=endt, length=500))
@test norm(y1[end]-y2[end])<0.15
@test LinearAlgebra.norm(y1[end]-y2[end])<0.15
end
println("ok.")

Expand Down
10 changes: 5 additions & 5 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using ODE, Test
using ODE, Test, LinearAlgebra

tols = [5e-2, 1e-2, 1e-2]

Expand Down Expand Up @@ -52,11 +52,11 @@ for solver in solvers
# -- = y ==> y = y0*e.^t
# dt
t,y=solver((t,y)->y, T(1), T[0:.001:1;])
@test maximum(abs.(y-e.^t)) < tol
@test maximum(abs.(y-.^t)) < tol
@test eltype(t)==T

t,y=solver((t,y)->y, T(1), T[1:-.001:0;])
@test maximum(abs.(y-e.^(t-1))) < tol
@test maximum(abs.(y .-.^ (t.-1))) < tol
@test eltype(t)==T

# dv dw
Expand All @@ -65,7 +65,7 @@ for solver in solvers
#
# y = [v, w]
t,y=solver((t,y)->[-y[2]; y[1]], T[1, 2], T[0:.001:2*pi;])
ys = hcat(y...).' # convert Vector{Vector{T}} to Matrix{T}
ys = hcat(y...)' # convert Vector{Vector{T}} to Matrix{T}
@test maximum(abs.(ys-[cos.(t)-2*sin.(t) 2*cos.(t)+sin.(t)])) < tol
@test eltype(t)==T
end
Expand All @@ -92,7 +92,7 @@ let
refsol = [0.2083340149701255e-07,
0.8333360770334713e-13,
0.9999999791665050] # reference solution at tspan[2]
@test norm(refsol-y[end], Inf) < 2e-10
@test LinearAlgebra.norm(refsol-y[end], Inf) < 2e-10
end
include("interface-tests.jl")
include("common.jl")
Expand Down

0 comments on commit 68227a1

Please sign in to comment.