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.

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

Fifth order EPIRK methods #417

Merged
merged 10 commits into from Jul 9, 2018
Merged

Fifth order EPIRK methods #417

merged 10 commits into from Jul 9, 2018

Conversation

MSeeker1340
Copy link
Contributor

Again, the integrators EPIRK5s3 and EXPRB53s3 comes from this paper. Unfortunately current EPIRK5s3 produces large test error and seems to have only first order convergence. This is most likely a result of a typo/bug and not because of the underlying framework (since the other fifth order method EXPRB53s3 works without problem).

# Compute the second group for U3
B = [zeros(R2) zeros(R2) zeros(R2) R2]
K = phiv_timestep([dt/2, 9dt/10], A, B; kwargs...)
U3 .+= 216/(25*dt^2) .* K[:, 1] + 8/dt^2 .* K[:, 2]
Copy link
Member

Choose a reason for hiding this comment

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

For my own future reference, 216 = 8*27, the 8 = 1/2^3 is the factor from the Krylov subspace algorithm because this calculates phi_i / t^i

Copy link
Member

Choose a reason for hiding this comment

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

Note that this line isn't fusing because of the scalar operations BTW. Just @. it. But optimizations can come later.

@ChrisRackauckas
Copy link
Member

I'm not sure there's a typo in your code. I double checked your coefficient calculations and it looked solid to me.

@MSeeker1340
Copy link
Contributor Author

I wrote a classical ExpRK implementation of EPIRK5s3 (using single step phiv instead of phiv_timestep) here: https://github.com/MSeeker1340/OrdinaryDiffEq.jl/commit/cb14d529f8f18039d34331941c8fe453ede9a9a9. For the same test problem in the test script, the two versions produce about the same error (~10^-4) whereas other ExpRK methods, be it classical or EPIRK, gives about 10^-7 error.

I have high suspicion that there is a mistake in the original formula. Since the coefficients are originally derived by symbolically solving a set of order conditions using Mathematica by the authors, I can try to substitute back the coefficients to check consistency.

Another angle I can tackle is to implement more EPIRK methods, like the original EPIRK5-P1/P2 integrators by Tokman. If all these don't have problems, then the likelihood that EPIRK5s3 is wrong will certainly increase.

@ChrisRackauckas
Copy link
Member

Please do EPIRK5-P1. That one seems to be the most efficient when it doesn't have order loss (i.e. not too stiff) so it would be a pretty essential algorithm to have. EPIRK5s3 is really beat out in any case by EXPRB53s3 anyways, so if it doesn't exist then it doesn't have much of a practical difference (though I do like the completionism to be able to verify benchmarks)

Copy link
Member

@ChrisRackauckas ChrisRackauckas left a comment

Choose a reason for hiding this comment

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

LGTM. Merge when tests pass, and depwarns can be fixed later.

@MSeeker1340
Copy link
Contributor Author

Note on the choice coefficients for EPIRK5-P1/P2: the $p$ coefficients aren't really free as they are cancelled out in the final form of the formula. The only free coefficient is $g_{22}$ in both cases. I made the simplest choice of $g_{22} = g_{32}$, so that no intermediate output is needed for the second column. I'm not sure if this is the optimal choice, and the paper doesn't have a discussion about this either.

@MSeeker1340 MSeeker1340 changed the title WIP: Fifth order EPIRK methods Fifth order EPIRK methods Jul 9, 2018
@coveralls
Copy link

coveralls commented Jul 9, 2018

Coverage Status

Coverage increased (+0.5%) to 74.809% when pulling 173fefc on MSeeker1340:expRK into c41ad05 on JuliaDiffEq:master.

@ChrisRackauckas ChrisRackauckas merged commit 6dbdf11 into SciML:master Jul 9, 2018
@ChrisRackauckas ChrisRackauckas mentioned this pull request Jul 15, 2018
@MSeeker1340 MSeeker1340 deleted the expRK branch July 17, 2018 01:17
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants