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

Conversation

sonVishal
Copy link

Hi,
I have tried to work out the component-wise relative and absolute tolerances for the non-stiff solvers based on 'Solving Ordinary Differential Equations I by Hairer and Wanner' and component-wise absolute tolerance combined with strictly scalar relative tolerance for the stiff solvers based on the MATLAB ode23s solver.
Since the interface checks does not handle vectors yet, the interface checks fail.
I have also changed the value of the error for the Robertson test since even Matlab's ode23s cannot produce a result with error strictly less than 2e-10. I got 2.0423543e-10 using the same options as used over here. Hence the change.
Please review.

Thanks and regards,
Vishal Sontakke,
Technical University of Munich.

…ing Ordinary Differential Equations I by Hairer and Wanner.

Component-wise AbsTol for Implicit Solvers similar to MATLAB ode23s.
tau = max(reltol*norm(x0, Inf), abstol)
d0 = norm(x0, Inf)/tau
tau = max(reltol.*abs(x0), abstol)
d0 = norm(x0,./tau Inf)
Copy link
Contributor

Choose a reason for hiding this comment

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

This looks wrong. Should probably be x0./tau, Inf.

Copy link
Author

Choose a reason for hiding this comment

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

Ah yes. I had fixed that in my code. Didn't update. Sorry. Will commit it again.

@acroy
Copy link
Contributor

acroy commented Mar 31, 2016

Thanks for the PR! You have accidentally deleted the .travis.yml file, which prevents us from getting Travis feedback. Please put it back.

@sonVishal
Copy link
Author

I added the file and also fixed the error as you pointed out.

# 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.

Now using only norm() since it is provided by the user and default value is provided.
@acroy
Copy link
Contributor

acroy commented Apr 1, 2016

Is there a fundamental reason why ode23s cannot have component-wise relative errors? I think one of the goals of this package is providing a uniform interface for different solvers (if possible).

Could we use a function for calculating the error estimate and facilitate mult. dispatch to distinguish between scalar and vector tolerances? That way we could keep the the old behavior for the scalar-like input and also avoid creating a tolerance vector for scalar abstol.

In case of scalar reltol "and" abstol, restore old behaviour.
@sonVishal
Copy link
Author

I couldn't find a fundamental reason behind such a restriction.

In theory reltol can also be specified component-wise.

I was able to implement component-wise reltol for ode23s similar to the non-stiff solvers.

I have tried to restore the old behaviour when reltol and abstol are both scalars.
Let me know if another function is really necessary for this.

In case of scalar reltol "and" abstol, restore old behaviour.

Also, passing "norm" to hinit so that the user defined norm can be used for the initial step size computation.
Should execute the old behaviour only if both abstol and reltol are scalars.
If either of them are vectors, perform component-wise computations.
@acroy
Copy link
Contributor

acroy commented Apr 1, 2016

Provided type inference works (reltol and abstol are keywords), separate functions and multiple dispatch allow you to avoid the if conditionals, which might have an impact on performance.

@sonVishal
Copy link
Author

Ah yes. I will try to update one with a function to get the error using multiple dispatch. 👍

@coveralls
Copy link

Coverage Status

Coverage decreased (-1.9%) to 91.768% when pulling d68c1a0 on sonVishal:master into c12fc36 on JuliaLang:master.

@ChrisRackauckas
Copy link
Member

What's the current status on this?

@pwl
Copy link

pwl commented Sep 2, 2016

We haven't planned for this in #49, but we should keep this in mind as a long term goal (after we merge #49). I think we should close this for now and open an issue instead.

@ChrisRackauckas
Copy link
Member

I don't think this should be closed due to an upcoming PR which doesn't include the functionality. It looks to me like this is complete, though I don't see a test. If @sonVishal adds some tests to show that it's working properly I'd vote to merge. How it relates with #49 is a separate issue.

@sonVishal
Copy link
Author

I will write a test case. I tested it myself but will include a test. Give me a week or so since I am currently very busy working on 2 different projects.

@sonVishal
Copy link
Author

I simply changed the relative and absolute tolerances for the ROBER test case to arrays and the test passed on my system.

Let me know if there is a need for another test case.

@coveralls
Copy link

coveralls commented Sep 18, 2016

Coverage Status

Changes Unknown when pulling f256c4d on sonVishal:master into * on JuliaDiffEq:master*.

@codecov-io
Copy link

Current coverage is 92.07% (diff: 80.76%)

No coverage report found for master at c12fc36.

Powered by Codecov. Last update c12fc36...f256c4d

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

8 participants