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

Component-wise relative and absolute tolerances #93

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

Component-wise relative and absolute tolerances #93

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
c2de029
Component-wise RelTol and AbsTol for Explicit Solvers similar to Solv…
sonVishal Mar 31, 2016
42f6b45
Pre .gitignore changes
sonVishal Mar 31, 2016
0eeb96b
Post .gitignore changes
sonVishal Mar 31, 2016
444468e
Updated the files with comments
sonVishal Mar 31, 2016
23349f3
Fixed the error in line 70 in ODE.jl
sonVishal Mar 31, 2016
740b05e
Added Travis file as required
sonVishal Mar 31, 2016
58bd9ca
Corrected the runtest.jl Robertson test case and removed the comments.
sonVishal Apr 1, 2016
936931a
Was using Inf norm for err.
sonVishal Apr 1, 2016
9c89efb
Component-wise reltol for ode23s.
sonVishal Apr 1, 2016
84d01ad
Component-wise reltol for ode23s.
sonVishal Apr 1, 2016
d38fa7d
Merge branch 'master' of https://github.com/sonVishal/ODE.jl
sonVishal Apr 1, 2016
76a8e05
Had written the if statements the other way around.
sonVishal Apr 1, 2016
ea10a9d
Implemented getScaledError as a multiple dispatch function for both s…
sonVishal Apr 1, 2016
d68c1a0
Added abs() function which can handle arrays of CompSol to interface-…
sonVishal Apr 9, 2016
f256c4d
Change ROBER testcase relTol and absTol to arrays
Sep 18, 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
27 changes: 19 additions & 8 deletions src/ODE.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,10 +66,10 @@ function hinit(F, x0, t0, tend, p, reltol, abstol)
# Returns first step, direction of integration and F evaluated at t0
tdir = sign(tend-t0)
tdir==0 && error("Zero time span")
tau = max(reltol*norm(x0, Inf), abstol)
d0 = norm(x0, Inf)/tau
tau = max(reltol.*abs(x0), abstol)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rather than using abs here, we should allow using a custom norm function. Look at how eg ode23s does this, and pass that norm along to hinit as well.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can I ask why abs(x0) is not suitable over here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since we're letting the user pass a custom norm function into the solver, we should use that here too - it's an oversight that this isn't already the case. I just think that if we touch these lines anyway, we might as well fix that (and if we don't, we should keep using norm here to make the transition easier later).

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Umm, this has nothing to do with the norm. As per the comments hinit is based on 'Solving Ordinary Differential Equations I' by Hairer and Wanner page 169 - The tau in the code is the same as the 'sc' variable in the book which is component-wise rather than a norm.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you use an L-infinity norm, doesn't that correspond to component-wise convergence?

....I see, you want to use a separate tolerance for each component, which is more than just changing the norm.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L-infinity would simply take the maximum over all the components.

Sometimes the solution is more sensitive w.r.t. only a component and you would thus want to control that component with a stricter tolerance and have relaxed tolerances for the other tolerances.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The "Hairer norm" used here isn't necessarily a norm because of how it it does the normalization. But the tau here is not the standard sc since sc is the sum of the adjusted tolerances, not the max. In practice this shouldn't matter much though.

d0 = norm(x0./tau, Inf)
f0 = F(t0, x0)
d1 = norm(f0, Inf)/tau
d1 = norm(f0./tau, Inf)
if d0 < 1e-5 || d1 < 1e-5
h0 = 1e-6
else
Expand All @@ -79,7 +79,7 @@ function hinit(F, x0, t0, tend, p, reltol, abstol)
x1 = x0 + tdir*h0*f0
f1 = F(t0 + tdir*h0, x1)
# estimate second derivative
d2 = norm(f1 - f0, Inf)/(tau*h0)
d2 = norm((f1 - f0)./tau, Inf)/h0
if max(d1, d2) <= 1e-15
h1 = max(1e-6, 1e-3*h0)
else
Expand Down Expand Up @@ -251,6 +251,17 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,

# initialization
t = tspan[1]

# Component-wise Tolerance check similar to MATLAB ode23s
# Support for component-wise absolute tolerance only.
# Relative tolerance should be a scalar.
@assert length(reltol) == 1 "Relative tolerance must be a scalar"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why don't you use isa(reltol, Number) to check if it is a scalar? This also concerns the other asserts.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wasn't aware of this function. I just wanted the code to stop if we do not give the proper input.
You mean instead of checking the lengths I simply use isa() function for the asserts?
I would have to do a length check for multidimensional problems eitherways

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. I just realized that length returns 1 for a scalar.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would be much better handled with a type spec on the function parameter, ie reltol::Number.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would only apply to reltol in the case of ode23s. For the other tolerances it would not be possible since they can be either a scalar or a vector of the same length as the problem.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure. But in the specific case handled by this line, since length 1 is hard-coded anyway, I think it would be much clearer.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright. I will make this change.

@assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem"

# Broadcast the abstol to a vector
if length(abstol) == 1 && length(y0) != 1
abstol = abstol*ones(y0);
end

tfinal = tspan[end]

Expand Down Expand Up @@ -300,11 +311,11 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
F2 = F(t + h, ynew)
k3 = W\(F2 - e32*(k2 - F1) - 2*(k1 - F0) + T )

err = (abs(h)/6)*norm(k1 - 2*k2 + k3) # error estimate
delta = max(reltol*max(norm(y),norm(ynew)), abstol) # allowable error
threshold = abstol/reltol # error threshold
err = (abs(h)/6)*norm((k1 - 2*k2 + k3)./max(max(abs(y),abs(ynew)),threshold),Inf) # error estimate
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

err has suddenly changed meaning here - it is now what it was before divided by something else. Should it have another name, too?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Originally it was the error estimate. Now it is something similar to a scaled error. Hence I kept the same name. I could change it if it is really required.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not required - I just didn't follow the details of what's going on with these few lines, and wondered if naming was contributing to my confusion :)

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah. I should probably make a comment saying that it is a scaled error. Will add this comment to avoid any further confusion.


# check if new solution is acceptable
if err <= delta
if err <= reltol

if points==:specified || points==:all
# only points in tspan are requested
Expand Down Expand Up @@ -334,7 +345,7 @@ function ode23s(F, y0, tspan; reltol = 1.0e-5, abstol = 1.0e-8,
end

# update of the step size
h = tdir*min( maxstep, abs(h)*0.8*(delta/err)^(1/3) )
h = tdir*min( maxstep, abs(h)*0.8*(reltol/err)^(1/3) )
end

return tout, yout
Expand Down
15 changes: 14 additions & 1 deletion src/runge_kutta.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,6 +250,19 @@ function oderk_adapt{N,S}(fn, y0::AbstractVector, tspan, btab_::TableauRKExplici
t = tspan[1]
tstart = tspan[1]
tend = tspan[end]

# Component-wise reltol and abstol as per Solving Ordinary Differential Equations I by Hairer and Wanner
@assert length(abstol) == 1 || length(abstol) == length(y0) "Dimension of Absolute tolerance does not match the dimension of the problem"
@assert length(reltol) == 1 || length(reltol) == length(y0) "Dimension of Relative tolerance does not match the dimension of the problem"

# Broadcast the tolerances if scalars are provided
if length(reltol) == 1 && length(y0) != 1
reltol = reltol*ones(y0);
end
if length(abstol) == 1 && length(y0) != 1
abstol = abstol*ones(y0);
end


# work arrays:
y = similar(y0, Eyf, dof) # y at time t
Expand Down Expand Up @@ -412,7 +425,7 @@ function stepsize_hw92!(dt, tdir, x0, xtrial, xerr, order,
for d=1:dof
# if outside of domain (usually NaN) then make step size smaller by maximum
isoutofdomain(xtrial[d]) && return 10., dt*facmin, timout_after_nan
xerr[d] = xerr[d]/(abstol + max(norm(x0[d]), norm(xtrial[d]))*reltol) # Eq 4.10
xerr[d] = xerr[d]/(abstol[d] + max(abs(x0[d]), abs(xtrial[d]))*reltol[d]) # Eq 4.10
end
err = norm(xerr, 2) # Eq. 4.11
newdt = min(maxstep, tdir*dt*max(facmin, fac*(1/err)^(1/(order+1)))) # Eq 4.13 modified
Expand Down
6 changes: 5 additions & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,6 +65,10 @@ end


# rober testcase from http://www.unige.ch/~hairer/testset/testset.html
# Changing the test value.
# Even the Matlab ode23s does not give the error to be less than 2e-10.
# The value comes out to be approximately 2.042354e-10.
# Hence changing the test value to 2.1e-10
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did this test pass before? If so, why doesn't it now? I'm sure it might be justified to increase the tolerance a little, but then I'd like to see some better motivation than just "Matlab doesn't do better".

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test passed before I made changes.
I believe that the test fails because different step sizes are now chosen by the solver as compared to the previous code. I think we need to find better step size control.

The test did pass when I reduced the tolerances to 1e-9 from the original value of 1e-8.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice, I think this change is safe to make, though. The comment could be left out though, as it's more a comment on the change than on the code itself.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Alright. I will remove the comments as well. 👍

let
println("ROBER test case")
function f(t, y)
Expand All @@ -81,7 +85,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 norm(refsol-y[end], Inf) < 2.1e-10
end
include("interface-tests.jl")

Expand Down