# Adaptivity in MIRK4 for BVPs

**URL:** <https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341>\
**Category:** Specific Domains\
**Tags:** diffeq\
**Created:** [October 1, 2019, 12:55am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341 "2019-10-01T00:55:30Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 1, 2019, 12:55am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/1 "2019-10-01T00:55:30Z")

</div>

I’m just messaging to ask if there is any plans to implement an adaptive mesh in the ODE-BVP solvers in BoundaryValueDiffEq.jl? I’m working on a problem where this is necessary, and for which the adaptive mesh in bvp4c in Matlab works quite well. Even the ability to specify a non-uniform mesh (rather than only specifying the dt value for a uniform mesh) would be very helpful.

Thanks

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 1, 2019, 7:01pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/2 "2019-10-01T19:01:40Z")

</div>

It’s planned… we just never got back to it. To teach @yingboma how adaptivity and interpolations of ODEs work, I had him do some ODE stuff and then he started optimizing our stiff ODE solvers to benchmark super well haha. We should return to that soon since it’s really not difficult and it’s just waiting to be done… (TBH the main reason is probably because no one on the dev team has very many BVPs to solve, so if you give a good reason we might whip adaptivity together 🙂)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 1, 2019, 7:37pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/4 "2019-10-01T19:37:24Z")

</div>

We have tons of them (almost 100?). This is specifically asking for an adaptive MIRK method for BVPs.

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [October 1, 2019, 7:38pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/5 "2019-10-01T19:38:01Z")

</div>

Isn’t that what approxfun is good at?

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 1, 2019, 7:38pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/6 "2019-10-01T19:38:05Z")

</div>

Haha, it’s good to hear the backstory. DifferentialEquations.jl is incredible - thanks for all your work. I thought of implementing adaptivity myself - as you say, conceptually it’s straightforward enough. However, a brief look through the source code for DifferentialEquations.jl has reminded me I might be more of a Julia user than a developer ☹

You asked for a planned use-case, so… 🙂

In chemical engineering, BVPs are arguably more common than IVPs. This is because we’re often interested in modelling spatially non-uniform physical systems, and in such cases you tend to get BCs in multiple locations simultaneously. For example, when designing distillation columns or absorbers, we’ll often specify some concentrations at the top and some at the bottom, and solve for the concentrations everywhere in between. Similarly, when dealing with reactive absorption of a gas into a liquid, we’ll have some boundary conditions at the surface of the liquid, and some from the bulk of the liquid. Secretly, I’d like to combine both into one complete model - right now such multi-scale models (simultaneously modelling a whole column and also of the local gas-liquid reactive dynamics) are only available in expensive, commercial solvers, where special care is taken to discretise everything just right and form carefully-structured matrices which are solved via Fortran. I suspect a julia implementation built off DifferentialEquations.jl could be conceptually much simpler and fast enough to be of practical use.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 1, 2019, 7:39pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/7 "2019-10-01T19:39:35Z")

</div>

> [@dpsanders](#):
>
> Isn’t that what approxfun is good at?

If it’s a smooth BVP, then yes ApproxFun.jl is a good option, and probably the best option.

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 1, 2019, 7:40pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/8 "2019-10-01T19:40:23Z")

</div>

That’s an interesting thought. Can approxfun handle stiff systems? I’ve found that Matlab’s bvp4c is particularly capable in that regard.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 1, 2019, 7:42pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/9 "2019-10-01T19:42:57Z")

</div>

> [@anodos](#):
>
> are only available in expensive, commercial solvers, where special care is taken to discretise everything just right and form carefully-structured matrices which are solved via Fortran.

Yeah, and that’s the part we cared to do. It has special forms for two-point BVPs so that it can utilize BandedMatrix types and then exploit the banded QR algorithm in that case.

> [@anodos](#):
>
> Can approxfun handle stiff systems? I’ve found that Matlab’s bvp4c is particularly capable in that regard.

Anything that isn’t shooting should do well there, since BVPs are inherently implicit. If your problem is ill-conditioned and smooth, spectral discretizations probably do better than MIRK methods.

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 1, 2019, 7:45pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/10 "2019-10-01T19:45:14Z")

</div>

Thanks for your help - I will give ApproxFun,jl a try.

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 1, 2019, 9:49pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/11 "2019-10-01T21:49:25Z")

</div>

Just to follow up on this - do you know if it is possible to define domains with breakpoints in ApproxFun.jl, as discussed here for Chebfun:

[https://www.chebfun.org/examples/ode-linear/Breakpoints.html](https://www.chebfun.org/examples/ode-linear/Breakpoints.html)

This would probably be ‘close enough’ to adaptive meshing for boundary layer problems. I’ve played around with the various ‘PiecewiseSegment’ options, but they don’t seem to play nice with the newton solver in ApproxFun.jl.  
Thanks

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [October 1, 2019, 10:40pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/12 "2019-10-01T22:40:51Z")

</div>

Cc @dlfivefifty

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [October 2, 2019, 5:56am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/13 "2019-10-02T05:56:41Z")

</div>

In theory it should work, though it’s combining two “proof of concept” pieces, so will require a number of bug fixes. Unfortunately I’m in the middle of completely revamping the underpinnings of ApproxFun so that it scales better (including to fully support hp methods) so I won’t have time to look at this, though any PR will be greatly appreciated.

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 2, 2019, 4:03pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/14 "2019-10-02T16:03:25Z")

</div>

Which two proof of concept pieces would need to be combined? I might have a look - I understand the theory, but don’t really understand Julia 🙂

---

<div class="post-metadata">

**Author:** ![morningglory05](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/morningglory05/32/43282_2.png) [@morningglory05](https://discourse.julialang.org/u/morningglory05)\
**Post date:** [February 1, 2023, 8:30pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/15 "2023-02-01T20:30:37Z")

</div>

Hello, I also would like to use a specified mesh rather than just defining a dt. Did this update already happen? If so, could you point me to documentation on it? Thanks!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 6, 2023, 12:40pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/16 "2023-02-06T12:40:06Z")

</div>

I’m not sure it takes a vector dt at this time.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 6, 2023, 1:43pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/17 "2023-02-06T13:43:38Z")

</div>

There are many different BVP solvers in BifurcationKit.jl. However, they are tailored for periodic orbits. Changing the boundary condition would not take too much effort. For example, you would need to change this [line](https://github.com/bifurcationkit/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitCollocation.jl#L445) for the collocation method.

You can specify the time mesh for all of them.

You can automatically adapt the time mesh for the orthogonal collocation algo during continuation. I would say, with a newton callback, you could perhaps adapt during newton iterations.

However, it is not ready for use as is for your use case.

---

<div class="post-metadata">

**Author:** ![anodos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/anodos/32/212778_2.png) [@anodos](https://discourse.julialang.org/u/anodos)\
**Post date:** [October 16, 2024, 2:33am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/18 "2024-10-16T02:33:23Z")

</div>

Just going to jump in and say it looks like ODEBVPs in Julia have undergone a massive upgrade! I assume this is related to,

> **[Google Summer of Code](https://summerofcode.withgoogle.com/archive/2023/projects/zLtIiHjs)**
>
> Google Summer of Code is a global program focused on bringing more developers into open source software development.

The library now has adaptive mesh refinement, interfaces with some solid Fortran libraries, as well as built-in support for multiple shooting methods (very, very handy for problems with challenging stability; I won’t show you the hand-written, 50-point multiple shooting method I once wrote for a widely-shared distillation column model.)

A big thanks for all this work!!! I’m in the process of putting together a 3rd year course on the modelling of chemical reactors, and this might just be the impetus to switch from Python to Julia.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [October 16, 2024, 5:53am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/19 "2024-10-16T05:53:25Z")

</div>

where did you see this?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2024, 10:41am UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/20 "2024-10-16T10:41:50Z")

</div>

Its in the Boundaryvaluediffeq part of the doc. We will do a write up soon. Indeed it has been a big push over the last year and there are tons of hidden improvements to mention

---

<div class="post-metadata">

**Author:** ![nk12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nk12/32/213131_2.png) [@nk12](https://discourse.julialang.org/u/nk12)\
**Post date:** [November 2, 2024, 12:18pm UTC](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341/21 "2024-11-02T12:18:50Z")

</div>

Hi!  
Do I get that correctly that adaptive non-uniform meshing is now available in BoundaryValueDiffEq?

[Next page](https://discourse.julialang.org/t/adaptivity-in-mirk4-for-bvps/29341.md?page=2)