Rkv76iia

[KeshavVenkatesh]

Checklist

• Added the algorithm for the RKV76IIA method
• Added a convergence test for the RKV76IIA method
• Creating this pull request to get feedback on the implementation

Additional Notes

• This particular convergence test fails during the order of convergence check.
• Info on failed tests:

1. Order between dt=1//8 and dt=1//16: -0.7854954880938649
   RKV76IIa Convergence Tests: Test Failed
   Expression: ≈(order, 7, atol = testTol)
   Evaluated: -0.7854954880938649 ≈ 7 (atol=0.3)

2. Order between dt=1//16 and dt=1//32: 1.5849625007211563
   RKV76IIa Convergence Tests: Test Failed at C:\Users\srira\Downloads\Ordinary Differential Equations\OrdinaryDiffEq.jl\lib\OrdinaryDiffEqVerner\test\rkv76iia_tests.jl:73
   Expression: ≈(order, 7, atol = testTol)
   Evaluated: 1.5849625007211563 ≈ 7 (atol=0.3)

3. Order between dt=1//32 and dt=1//64: 0.967333811079678
   RKV76IIa Convergence Tests: Test Failed at C:\Users\srira\Downloads\Ordinary Differential Equations\OrdinaryDiffEq.jl\lib\OrdinaryDiffEqVerner\test\rkv76iia_tests.jl:73
   Expression: ≈(order, 7, atol = testTol)
   Evaluated: 0.967333811079678 ≈ 7 (atol=0.3)

[ChrisRackauckas]

Show plots of the convergence tests, just plot(sim) using the plot recipe

[KeshavVenkatesh]

,

Show plots of the convergence tests, just plot(sim) using the plot recipe

Could you please let me know if this is the plot you are looking for?

[KeshavVenkatesh]

I will also take a look at the RKV76IIA algorithm's implementation to find out why the order of convergence is failing at smaller values of dt. If you have any suggestions as to what I should look for specifically, please let me know.

[oscardssmith]

Your plot is using too small dt for Float64 tolerance. you can see that it levels off at y=~1e-14 which is roughly floating point error. You should probably switch to using BigFloat for the test.

[ChrisRackauckas]

Yes it's just saturating because that is as accurate as Float64 can get. Doing this test in big float would let it keep going

[KeshavVenkatesh]

I tried to convert the errors obtained from the convergence test to BigFloat, but that did not help with the convergence. This is the output for the code:

Out-of-place solution at t=1: 0.36787944117147253
Expected value: 0.36787944117144233

Testing order 7:
dt = 1//4, error = 4.087773236764922e-12
dt = 1//8, error = 3.798205619585314e-14
dt = 1//16, error = 6.546007605532576e-14
dt = 1//32, error = 2.1828311187557115e-14
dt = 1//64, error = 1.1170170151155611e-14
Order between dt=1//4 and dt=1//8: 6.749853347562831925005443312238533940673969705963336027706295162463012652846084
Order between dt=1//8 and dt=1//16: -0.7852972692785772597867943127951613256220467034287479973075014653126841550382554
Order between dt=1//16 and dt=1//32: 1.584414762352023612766097130426965010072873882633238027110177664600350876210189
Order between dt=1//32 and dt=1//64: 0.9665493541991219980382499685074483282026106254799138626449239186346385050363408
Test Summary: Total Time
RKV76IIa Convergence Tests 0 2.8s

I also looked at some tests defined for OrdinaryDiffEqLowStorageRK, and I tried to write similar tests for the RKV76IIA method, but the code is generating this compilation error:

RKV76IIa: Error During Test at OrdinaryDiffEq.jl/lib/OrdinaryDiffEqVerner/test/convergence_tests.jl:64
Got exception outside of a @test
MethodError: no method matching OrdinaryDiffEqVerner.RKV76IIaTableau(::Type{BigFloat}, ::Type{BigFloat})

As far as I can see, it is not easy to convert the RKV76IIA algorithm to support BigFloat without major refactoring.

Could you please provide specific guidance on how I can resolve the BigFloat issue and get the tests working? I am also attaching the code I have written as reference (convergence_test.txt).

[ChrisRackauckas]

This dispatch will only get the numbers in Float64, which is why it cuts off. For BigFloats dispatch, write it like this:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl#L892-L902

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl#L458-L465

The reason is that machine floats are only 64-bit, so 0.01710144927536231884057971014492753623188 will automatically truncate to the Float64. You need to input it as a string and tell it to parse the string to a BigInt/BigFloat in order to force it to keep the full precision. That requires a separate dispatch because it's slower, so this one is hit for CompiledFloats (i.e. Float32, Float64) while the other will be for arbitrary other number types, like BigFloat.

[ChrisRackauckas]

And leave the dispatch open there https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl#L438

[KeshavVenkatesh]

Thanks for your suggestions. I incorporated them to the tableau file and repeated the convergence tests.

I converted the errors obtained from the algorithm to BigFloat, but that did not help with the convergence. I tested the errors and convergences for the Vern7 and RK67IIA algorithms. As you can see, both the algorithms show similar error patterns for smaller values of dt.

Algorithm = RKV76IIa
Expected value: 0.36787944117144233
Testing order 7:
dt = 0.5, error = 5.1965687308808128364034928381443023681640625e-10, sol = 0.36787944065178546
dt = 0.25, error = 4.189981694935340783558785915374755859375e-12, sol = 0.36787944116725235
dt = 0.125, error = 6.0784710598227320588193833827972412109375e-14, sol = 0.3678794411715031
dt = 0.0625, error = 4.6906922790412863832898437976837158203125e-14, sol = 0.36787944117148924
dt = 0.03125, error = 2.4258373088059670408256351947784423828125e-14, sol = 0.3678794411714666
dt = 0.015625, error = 1.1934897514720432809554040431976318359375e-14, sol = 0.36787944117145427
Order between dt=0.5 and dt=0.25: 6.954471581718476389439298600217739628517308528761066969346766453557227185836897
Order between dt=0.25 and dt=0.125: 6.107091647999272164571434293710687406083044961407325939988096760821410796558153
Order between dt=0.125 and dt=0.0625: 0.3739076233189891030338734217531834679697423885325258587722005348902813160993851
Order between dt=0.0625 and dt=0.03125: 0.9513180616689483344043098717825337917119788813568864435639971719799842337317316
Order between dt=0.03125 and dt=0.015625: 1.023296619911138113569409709235687090822943924735172497034207897329299288988322

Algorithm = Vern7
Expected value: 0.36787944117144233
Testing order 7:
dt = 0.5, error = 9.92350590589552439269027672708034515380859375e-10, sol = 0.3678794421637929
dt = 0.25, error = 8.904932347064686837256886065006256103515625e-12, sol = 0.36787944118034727
dt = 0.125, error = 6.8833827526759705506265163421630859375e-14, sol = 0.36787944117151117
dt = 0.0625, error = 4.44089209850062616169452667236328125e-16, sol = 0.3678794411714428
dt = 0.03125, error = 5.5511151231257827021181583404541015625e-17, sol = 0.3678794411714424
dt = 0.015625, error = 5.5511151231257827021181583404541015625e-17, sol = 0.3678794411714423
Order between dt=0.5 and dt=0.25: 6.80010144374439438425040709903372447327888923308563360318083561540784168937752
Order between dt=0.25 and dt=0.125: 7.015343106941082170173433681622139920488953302324273065461036730872138749459362
Order between dt=0.125 and dt=0.0625: 7.276124405274237800112088335184954886350303105758803571730235069808926118048688
Order between dt=0.0625 and dt=0.03125: 3.000000000000000100370320137762229410836612506849573639125476871813688904806905
Order between dt=0.03125 and dt=0.015625: 0.0

Please let me know if I should investigate further as to why the RKV76IIA method is not converging.

[oscardssmith]

try with dt=1,2,4? Those should be big enough (based on the 0.5 and 0.25 values) to see if it is converging as expected.

[KeshavVenkatesh]

I tried using larger values (from 4 to 1/64), and here is the output that was generated:

Out-of-place solution at t=12.0: 6.144212353057301e-6
Expected value: 6.14421235332821e-6

Algorithm = RKV76IIa
Testing order 7 on [0, 12.0]:
dt = 4.0, error = 0.0013678068737745335875356333769914396469857820193283259868621826171875, sol(t=12.0) = -0.0013616626614212054
dt = 2.0, error = 2.94764395085484929258379249716881531639955937862396240234375e-08, sol(t=12.0) = 6.114735913819661e-6
dt = 1.0, error = 2.25747100651151901516688891291551044560037553310394287109375e-11, sol(t=12.0) = 6.144189778618145e-6
dt = 0.5, error = 1.04139788765643566403884534565804642625153064727783203125e-13, sol(t=12.0) = 6.144212249188421e-6
dt = 0.25, error = 8.3637049651426320640013045704108662903308868408203125e-16, sol(t=12.0) = 6.144212352491839e-6
dt = 0.125, error = 1.26182498152473121511008002926246263086795806884765625e-17, sol(t=12.0) = 6.144212353340828e-6
dt = 0.0625, error = 9.1352503361376291568518581698299385607242584228515625e-18, sol(t=12.0) = 6.144212353337345e-6
dt = 0.03125, error = 4.812841206298934526586208448861725628376007080078125e-18, sol(t=12.0) = 6.144212353333023e-6
dt = 0.015625, error = 2.427596426830824771769812286947853863239288330078125e-18, sol(t=12.0) = 6.144212353330637e-6

Order between dt = 4.0 and dt = 2.0: 15.50194274880805265698956720433776189758202794918343955342964318487648014191817
Order between dt = 2.0 and dt = 1.0: 10.35063909611411833028441194190386330160317414268724292473790224867003907536299
Order between dt = 1.0 and dt = 0.5: 7.760042263043809601916854423634858922374767836952699334292038210859350815195202
Order between dt = 0.5 and dt = 0.25: 6.960163499315261686143930957572082034960722078318739852404972795250671524110996
Order between dt = 0.25 and dt = 0.125: 6.050558447189970656003510831003772066007298706124645614383212729595953282820191
Order between dt = 0.125 and dt = 0.0625: 0.4659956494310546374988981040928102789602508885644379961104698796868095862405011
Order between dt = 0.0625 and dt = 0.03125: 0.9245554398582751924878797230142207918487469110046569766592370648854372865303651
Order between dt = 0.03125 and dt = 0.015625: 0.9873602219273955398789573300533855976925713386725199556123231753501798792685471

The orders are still very off for values larger than 1 and values smaller than 0.25. I will try to look at the RKV76IIA method tomorrow to see if there is anything wrong with it.

[ChrisRackauckas]

Share the current plot

[KeshavVenkatesh]

Share the current plot

[ChrisRackauckas]

One of the coefficients is wrong in the 20th digit or cutoff

[ChrisRackauckas]

Sum the b's and subtract 1. That should be zero. That should eyeball accuracy of one part. Then do the same with the rows of a to c

[KeshavVenkatesh]

One of the coefficients is wrong in the 20th digit or cutoff

Thank you

[KeshavVenkatesh]

I have made the following changes to my code:

1. I cleaned up all of the code in the testing.
2. I compared the RKV76IIa method with the Vern7 method for many different differential equations. The two had very similar outputs for the integration, which is good.
3. I checked the tableau file, and one of the values was off. Thanks for catching that; I have fixed it.
4. I checked my code again and made sure that all of the values were of type BigFloat.

Now, the order of the method is very close to 7 for small values. I believe this is due to some of the values not being BigFloats. I have attached all of the output for my code (testing the method against the Vern7 method and for different equations) below. You can find the orders and the convergence plot at the very bottom.

Sorry for all the delays in making this. Can you please look at the results and let me know how we should proceed next?

u(t) = sin(t)

Comparing dt usage across algorithms

In Place

Expected Solution - 0.8415

Computed Solution:
Algorithm: RKV76IIa
Solution at t=1.0: 0.8414709848078934
Final dt: 0.05952538518208472

Algorithm: Vern7
Solution at t=1.0: 0.8414709848078968
Final dt: 0.04015799693304212

Algorithm: SSPRK22
Solution at t=1.0: 0.7701511529340699
Final dt: 1.0

Comparing dt usage across algorithms

Out of Place

Expected Solution - 0.8415

Algorithm: RKV76IIa
Solution at t=1.0: [0.8414698511088894]
Final dt: 3.3370486541839384e-6

Algorithm: Vern7
Solution at t=1.0: [0.8414709848078968]
Final dt: 0.05901945378963047

Algorithm: SSPRK22
Solution at t=1.0: [0.7701511529340699]
Final dt: 1.0

du/dt = sin(10t) * u * (1 - u)

Comparing dt usage across algorithms

Expected Solution - ≈ 0.11899214697726988

Algorithm: RKV76IIa
Solution at t=1000.0: 0.1189921469744795
Final dt: 0.007386901132917956

Algorithm: Vern7
Solution at t=1000.0: 0.11899214702210924
Final dt: 0.006452126193153163

Algorithm: SSPRK22
Solution at t=1000.0: 1.0743872002703586e-126
Final dt: 1.0

u(t) = e ^ (-t)

Expected Solution 1.603810890548637926589339246984480554458308889336511692723506695181194395879304e-28

Comparing dt usage across algorithms

Algorithm: RKV76IIa
Solution at t=64.0: 1.6038108905494776e-28
Final dt: 0.02351369207259779

Algorithm: SSPRK22
Solution at t=64.0: 5.421010862427522e-20
Final dt: 1.0

Calling Vern7 shows this warning that I have to debug a little further
Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase ~/.julia/packages/SciMLBase/YE7xF/src/integrator_interface.jl:623

Algorithm: Vern7
Solution at t=1.2333480035129211e-18: 1.0
Final dt: 1.2332883520467626e-24

Testing convergence of RKV76IIa method

Computed Solution 1.6038108905653896e-28 for dt = 0.125, error = 1.675175597517680603479580639042602251693553433738067435466019026323254432543159e-39

Computed Solution 1.6038108893818258e-28 for dt = 0.25, error = 1.166812056131833020525879283060969634738675287642719275816700899152061110466407e-37

Computed Solution 1.6038107455651205e-28 for dt = 0.5, error = 1.449835173477254493009381262270676841176136295643344424930958373579210399672479e-35

Computed Solution 1.6037794634513104e-28 for dt = 1.0, error = 3.142709732742637762440293841727026139070200547295464662844181496811395567146923e-33

Computed Solution 1.563199573275837e-28 for dt = 2.0, error = 4.061131727280089672348812672221127635017627354724964466103919926683040411247637e-30

Computed Solution 5.188448212746488e-16 for dt = 4.0, error = 5.188448212744884511450189060011182528221470665612852283706204362617532505243073e-16

Computed Solution -5.008043695718157e7 for dt = 6.0, error = 5.00804369571815729141235351562500001603810890548637852976087034142335380998364e+07

Computed Solution 5.244491884184406e16 for dt = 8.0, error = 5.244491884184406399999999999999999999999999983961891094513621470239129658576616e+16

Order between dt=8.0 and dt=6.0: 29.96390870717052736845441586429084828076659924816814998210475561741697342794351

Order between dt=6.0 and dt=4.0: 76.3532902178153814052005827422066053641746152358764842839834119391110010053967

Order between dt=4.0 and dt=2.0: 46.8604146195877054703911030612798320678145912245210969142445524847132646164008

Order between dt=2.0 and dt=1.0: 10.33565708013339132898530903793513389399477045443897777018095630520321666989135

Order between dt=1.0 and dt=0.5: 7.75997632295133487103881686849804582617005129053679816493306817559997861727266

Order between dt=0.5 and dt=0.25: 6.957172886421031833799791174910677866045606883513609521182381582081417437125622

Order between dt=0.25 and dt=0.125: 6.122116056280128739688797786481122825680486891675880130192327941317864431186419

[ChrisRackauckas]

Okay yes that looks correct now. So fix the dt's for the convergence test into the range of acceptable values and it should be pass and that would be right.

[ChrisRackauckas]

Are these to sufficient accuracy? Is this exactly what the paper says? We can just compute them to match for example a31+a32 which might be more accurate.

[KeshavVenkatesh]

I checked the links referenced in this issue and added up each of the a values. It looks like a lot of the digits summed up to 10, so the c values didn't have very many digits.

[ChrisRackauckas]

This doesn't need to be separate. Integrate it into the other convergence tests, just make sure dt is set in the valid range.

[KeshavVenkatesh]

I wasn't able to find any convergence tests for the other algorithms in the OrdinaryDiffEqVerner directory. So, I integrated those into my RKV76IIa convergence test code.

I created a new file called 'ode_verner_tests.jl' in my pull request, which has the convergence tests for all the algorithms. However, I noticed that the Vern6 algorithm was failing. The overall convergence order was estimated to be 7.23868839418773, while the order was supposed to be 6. I have commented out the Vern6 algorithm in my code, but let me know if you want me to use a custom set of dts for this algorithm.

[KeshavVenkatesh]

If you have any other problems that I can work on, please let me know. I can work on a long-term project and simultaneously help you with smaller tasks as well.

[KeshavVenkatesh]

Okay yes that looks correct now. So fix the dt's for the convergence test into the range of acceptable values and it should be pass and that would be right.

I did this; here is the new plot:

[ChrisRackauckas]

why is this a new file instead of editing the previous one?

[KeshavVenkatesh]

The old file was just for testing the RKV76IIa algorithm out. I have removed that file from the PR now that I have written a more formal convergence test code.

[ChrisRackauckas]

The verner should probably be split off from

OrdinaryDiffEq.jl/test/regression/ode_unrolled_comparison_tests.jl

Line 104 in 2e0531a

tabalg = ExplicitRK(tableau = constructVernerEfficient6(BigFloat))

so the tests are moved to the new subrepo format.

[ChrisRackauckas]

Tests fail because they need SSPRK things... why are those there at all? Verner methods are not SSPRK methods so they shouldn't be testing any of that really

[KeshavVenkatesh]

Tests fail because they need SSPRK things... why are those there at all? Verner methods are not SSPRK methods so they shouldn't be testing any of that really

The old file (rkv76iia_tests.jl) was comparing the RKV76IIa algorithm to an algorithm in the OrdinaryDiffEqSSPRK directory. I have removed the rkv76iia_tests.jl file now.

[KeshavVenkatesh]

I took the code in the OrdinaryDiffEq.jl/test/regression/ode_unrolled_comparison_tests.jl file and moved it to the convergence test code (ode_verner_tests,jl). Could you please let me know if the convergence tests cover what you were looking for?

[ChrisRackauckas]

You need ExplicitRK to be a test dependency if you're going to use it here. That's going to be a CI failure.

[ChrisRackauckas]

Yes that test change looks good, though the Project.toml needs to be updated so there's a test dependency for the ExplicitRK method, I think it's just OrdinaryDiffEqExplicitRK?

[KeshavVenkatesh]

Yes that test change looks good, though the Project.toml needs to be updated so there's a test dependency for the ExplicitRK method, I think it's just OrdinaryDiffEqExplicitRK?

I updated Project.toml with the necessary changes, but I am still seeing the CI error. I am not sure if it is related to OrdinaryDiffEqVerner or not. Could you please let me know if I am missing anything?

[KeshavVenkatesh]

Can you please let me know if I need to fix anything in my code to make the error disappear? Thank you!

[ChrisRackauckas]

this isn't the right tableau, so I wouldn't expect it to exactly match yours.

[ChrisRackauckas]

should probably just remove the other solve and this check since it won't tableau match

[ChrisRackauckas]

Thanks!

[ChrisRackauckas]

Oh wait, just realized I forgot to remove those bits. Can you make a follow up to remove those extra tableaus?