5 ms·

A lazy way to solve differential equations

204 points by cosmic_quanta 4y ago

cosmic_quanta 4y ago
This post highlights why laziness can sometimes be a good thing. We talk a lot about space leaks and strange performance characteristics, but this post would not have been possible without infinite streams. Pretty cool stuff.
- freemint 4y ago
  This post is very possible without infinite streams. A for loop or a while loop would do the trick just as well.
  - cosmic_quanta 4y ago
    Of course, I'm not pretending the opposite. In practical terms, with finite precision, there no need for infinite streams. There are other comments in this thread about alternate approaches. Ignoring issues of convergence etc., the approach in the post can technically reach a level of precision approaching infinity. At the end of the day, this is just a neat demonstration of something which is made more elegant by infinite streams. Edit: I see the wording in my parent comment wasn't great. I meant: > ... but (the elegance of) this post would not have been possible without infinite streams. Pretty cool stuff.
- lupire 4y ago
johnbcoughlin 4y ago
The post describes essentially a spectral method for solving ODEs in the Taylor (monomial) basis. Haskell definitely makes it nice to play around with this, but the monomial basis is terribly ill conditioned and not a good choice for this application. The Chebyshev basis is much better. Anyone who finds this neat, I recommend checking out the chebfun family of packages. They use the same idea but with much more rigorous numerics.
- jamessb 4y ago
  > I recommend checking out the chebfun family of packages. There's also an introductory textbook about ODEs ("Exploring ODEs") that uses chebfun [2] throughout. As well as the original MATLAB implementation, there are now implementations in Julia, Python, and some other languages [3]. [1]: https://people.maths.ox.ac.uk/trefethen/ExplODE/ https://people.maths.ox.ac.uk/trefethen/ExplODE/ [2]: https://www.chebfun.org/ https://www.chebfun.org/ [3]: https://www.chebfun.org/about/projects.html https://www.chebfun.org/about/projects.html
- ChrisRackauckas 4y ago
  > The Chebyshev basis is much better The Chebyshev basis in time is still not great. This was used in Chebfun before but there's a note that for time to just use standard ODE solvers because it did not good performance. Collocation really only makes sense on BVPs.
  - mturmon 4y ago
    This comment and GP comment open up a helpful learning experience for a casual/occasional user like me. This distinction between initial value problems (IVP, one endpoint known, hence a “time” differential equation) and BVPs (both endpoints known, not “time”) also appears in this introduction to the Chebfun solvers: https://www.chebfun.org/docs/guide/guide07.html https://www.chebfun.org/docs/guide/guide07.html (last paragraph of sec 7.1) Following the breadcrumbs to the IVP treatment in chapter 10 https://www.chebfun.org/docs/guide/guide10.html https://www.chebfun.org/docs/guide/guide10.html (Sec 10.2) we read: A major change was introduced in Version 5.1 in how initial-value (and final-value) problems are solved. Before, Chebfun used the same global spectral representations as for BVPs. This usually works fine for linear problems, but for nonlinear ones, it is inferior to the method of time-stepping by Runge-Kutta or Adams formulas. In Chebfun Version 5.1, we accordingly switched to solving IVPs numerically by ode113 (by default), converting the resulting output to a Chebfun representation. Thanks again for the helpful comments.
- naillo 4y ago
  Anyone here (or the OP) know how you'd use chebyshev in a multivariate scenario? The wikipedia article only defines them for the 1d case.
lupire 4y ago
The OP paper is as simple as the blog post, laying out the conductive form for expressing the mathematical definitions of the college-level math concepts. So you can try porting the other sections to executable code also, as exercises. http://www.matematicas.unam.mx/favio/publ/cocalculus.pdf http://www.matematicas.unam.mx/favio/publ/cocalculus.pdf A top hit for more info on this topic, is this 2015 HN comment thread: https://news.ycombinator.com/item?id=9617083 https://news.ycombinator.com/item?id=9617083
superb-owl 4y ago
Nit: the FTC example has a typo - the integral should end at `a+x`, not `a`
- Blammar 4y ago
  Yes! For a while, I thought I had lost my mind, my decrepitude had finally arrived, and I was no longer remembering my teenage calculus.
freemint 4y ago
> analytic everywhere and when this hypothesis doesn’t hold, the method will just silently fail and return garbage such as infinity or NaN for the. This method will produce false results and not fail if the right hand side is not analytic but infinitely differentiable. It will also silently fail depending on the implementation of the automatic differentiation in the presence of integrated discontinuities anywhere in the function.
- whatshisface 4y ago
  I could be misunderstanding but it may also produce a series that converges, but poorly. When using Taylor series to solve differential equations, there's a final step of "try and recognize the function," so that you can have sin(x) instead of a power series that will take a billion terms to reasonably compute sin(100).
ChrisRackauckas 4y ago
In Julia, this is implemented with the number type of TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl and captured in the integration package TaylorIntegration.jl https://github.com/PerezHz/TaylorIntegration.jl https://github.com/PerezHz/TaylorIntegration.jl which is then available in the DifferentialEquations.jl common interface: https://docs.sciml.ai/dev/modules/DiffEqDocs/solvers/ode_solve/#TaylorIntegration.jl https://docs.sciml.ai/dev/modules/DiffEqDocs/solvers/ode_sol... The implementation described here is missing a lot of the optimizations that TaylorIntegration.jl employs (it makes the method no longer simple), but after doing some of the extra stuff you can find that Taylor integration methods can be good if you need really high accuracy. Some examples of this are seen in some of the SciMLBenchmarks, particularly https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_e... and https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.html https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_bos... . One of the nice things about Taylor integration methods is that they can be used to build interval solutions of ODEs which guarantee the solution. That's described in some detail in this video: https://www.youtube.com/watch?v=o1h7BUW04NI https://www.youtube.com/watch?v=o1h7BUW04NI .
- helmholtz 4y ago
  Hahaha, my first thought as I clicked the comments was "I wonder if Chris already has this in DifferentialEquations.jl" and here we are.
- joshjob42 4y ago
  I knew I'd seen something like this in Julia. Fantastic to see it's already implemented!
freemint 4y ago
I just noticed: > For example, do you want to approximate some complicated integral? Just use the Stream x that we previously defined: > ghci> erf = 0 :> exp (-x^2) > ghci> take 10 (toList erf) > [0.0,1.0,-0.0,-2.0,0.0,12.0,0.0,-120.0,0.0,1680.0] This is not even close to the tailor series of the erf function. This to me looks more like x' = exp(-x^2) which is entirely unrelated to the integral.
- deleted 4y ago
  [deleted]
- texaslonghorn5 4y ago
  The two are related because erf(x) is the integral of exp(-t^2) from 0 to x. So this is correct if you interpret coefficients in the basis {x^n/n!} and drop the scaling factor of 2/sqrt(pi). The Taylor series of exp has a 1 in the zeroth coefficient, so integrating from 0 to x we get a 1 in the 1st coefficient of erf. exp(-x^2) = 1 - x^2 + x^4/2 - x^6/6 + ... = 1/0! - 2x^2/2! + 12x^4/4! - 120x^6/6! + ... erf(x) = x/1! - 2x^3/3! + 12x^5/5! - 120x^7/7! + ...
snake_plissken 4y ago
Can someone better explain what the author means early on when they describe the Fundamental Theorem of Calculus as f(a + x) = f(a) + f'(t)dt|x -> a I thought that the integration yields the area under the curve from x to a? So how does the that relate to the output of f(a + x) which is not an area quantity?
- cat_man 4y ago
  I think it's a typo. \int_a^x f'(t)dt should be f(x) - f(a), and that would be the fundamental theorem of calculus (or some corollary of it). It looks like either f(a+x) should be f(x) or the upper limit x should be x+a. The rest of the post deals with Taylor series, which are often written with f(a+x) for the function being expanded, so that might have mixed the author up when writing that. They don't directly use that expression elsewhere as far as I can tell, so it doesn't affect the rest of the post. If it's f(x) = f(a) + \int_a^x f'(t)dt, then you can think of f(x) - f(a) as the area under the curve of f'(t) on the interval from a to x. For example, if f(x) = x, f'(x) = 1. The area of a function with height of 1 from a to x would be x - a, which is the same as f(x) - f(a) = x - a.
- deleted 4y ago
  [deleted]
gaze 4y ago
This seems a bit icky as you move away from the origin. There’s no promises for how fast the series converges and how well it does yea far out and so on and so forth. This reminds me a bit of harmonic balance which is a great, well behaved spectral technique

204 points by cosmic_quanta 4y ago

This post highlights why laziness can sometimes be a good thing. We talk a lot about space leaks and strange performance characteristics, but this post would not have been possible without infinite streams. Pretty cool stuff.

freemint 4y ago
This post is very possible without infinite streams. A for loop or a while loop would do the trick just as well.
- cosmic_quanta 4y ago
  Of course, I'm not pretending the opposite. In practical terms, with finite precision, there no need for infinite streams. There are other comments in this thread about alternate approaches. Ignoring issues of convergence etc., the approach in the post can technically reach a level of precision approaching infinity. At the end of the day, this is just a neat demonstration of something which is made more elegant by infinite streams. Edit: I see the wording in my parent comment wasn't great. I meant: > ... but (the elegance of) this post would not have been possible without infinite streams. Pretty cool stuff.
lupire 4y ago

johnbcoughlin 4y ago

The post describes essentially a spectral method for solving ODEs in the Taylor (monomial) basis. Haskell definitely makes it nice to play around with this, but the monomial basis is terribly ill conditioned and not a good choice for this application. The Chebyshev basis is much better. Anyone who finds this neat, I recommend checking out the chebfun family of packages. They use the same idea but with much more rigorous numerics.

jamessb 4y ago
> I recommend checking out the chebfun family of packages. There's also an introductory textbook about ODEs ("Exploring ODEs") that uses chebfun [2] throughout. As well as the original MATLAB implementation, there are now implementations in Julia, Python, and some other languages [3]. [1]: https://people.maths.ox.ac.uk/trefethen/ExplODE/ https://people.maths.ox.ac.uk/trefethen/ExplODE/ [2]: https://www.chebfun.org/ https://www.chebfun.org/ [3]: https://www.chebfun.org/about/projects.html https://www.chebfun.org/about/projects.html
ChrisRackauckas 4y ago
> The Chebyshev basis is much better The Chebyshev basis in time is still not great. This was used in Chebfun before but there's a note that for time to just use standard ODE solvers because it did not good performance. Collocation really only makes sense on BVPs.
- mturmon 4y ago
  This comment and GP comment open up a helpful learning experience for a casual/occasional user like me. This distinction between initial value problems (IVP, one endpoint known, hence a “time” differential equation) and BVPs (both endpoints known, not “time”) also appears in this introduction to the Chebfun solvers: https://www.chebfun.org/docs/guide/guide07.html https://www.chebfun.org/docs/guide/guide07.html (last paragraph of sec 7.1) Following the breadcrumbs to the IVP treatment in chapter 10 https://www.chebfun.org/docs/guide/guide10.html https://www.chebfun.org/docs/guide/guide10.html (Sec 10.2) we read: A major change was introduced in Version 5.1 in how initial-value (and final-value) problems are solved. Before, Chebfun used the same global spectral representations as for BVPs. This usually works fine for linear problems, but for nonlinear ones, it is inferior to the method of time-stepping by Runge-Kutta or Adams formulas. In Chebfun Version 5.1, we accordingly switched to solving IVPs numerically by ode113 (by default), converting the resulting output to a Chebfun representation. Thanks again for the helpful comments.
naillo 4y ago
Anyone here (or the OP) know how you'd use chebyshev in a multivariate scenario? The wikipedia article only defines them for the 1d case.

lupire 4y ago

The OP paper is as simple as the blog post, laying out the conductive form for expressing the mathematical definitions of the college-level math concepts. So you can try porting the other sections to executable code also, as exercises. http://www.matematicas.unam.mx/favio/publ/cocalculus.pdf http://www.matematicas.unam.mx/favio/publ/cocalculus.pdf A top hit for more info on this topic, is this 2015 HN comment thread: https://news.ycombinator.com/item?id=9617083 https://news.ycombinator.com/item?id=9617083

superb-owl 4y ago

Nit: the FTC example has a typo - the integral should end at `a+x`, not `a`

Blammar 4y ago
Yes! For a while, I thought I had lost my mind, my decrepitude had finally arrived, and I was no longer remembering my teenage calculus.

freemint 4y ago

> analytic everywhere and when this hypothesis doesn’t hold, the method will just silently fail and return garbage such as infinity or NaN for the. This method will produce false results and not fail if the right hand side is not analytic but infinitely differentiable. It will also silently fail depending on the implementation of the automatic differentiation in the presence of integrated discontinuities anywhere in the function.

whatshisface 4y ago
I could be misunderstanding but it may also produce a series that converges, but poorly. When using Taylor series to solve differential equations, there's a final step of "try and recognize the function," so that you can have sin(x) instead of a power series that will take a billion terms to reasonably compute sin(100).

ChrisRackauckas 4y ago

In Julia, this is implemented with the number type of TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl and captured in the integration package TaylorIntegration.jl https://github.com/PerezHz/TaylorIntegration.jl https://github.com/PerezHz/TaylorIntegration.jl which is then available in the DifferentialEquations.jl common interface: https://docs.sciml.ai/dev/modules/DiffEqDocs/solvers/ode_solve/#TaylorIntegration.jl https://docs.sciml.ai/dev/modules/DiffEqDocs/solvers/ode_sol... The implementation described here is missing a lot of the optimizations that TaylorIntegration.jl employs (it makes the method no longer simple), but after doing some of the extra stuff you can find that Taylor integration methods can be good if you need really high accuracy. Some examples of this are seen in some of the SciMLBenchmarks, particularly https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_e... and https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.html https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_bos... . One of the nice things about Taylor integration methods is that they can be used to build interval solutions of ODEs which guarantee the solution. That's described in some detail in this video: https://www.youtube.com/watch?v=o1h7BUW04NI https://www.youtube.com/watch?v=o1h7BUW04NI .

helmholtz 4y ago
Hahaha, my first thought as I clicked the comments was "I wonder if Chris already has this in DifferentialEquations.jl" and here we are.
joshjob42 4y ago
I knew I'd seen something like this in Julia. Fantastic to see it's already implemented!

freemint 4y ago

I just noticed: > For example, do you want to approximate some complicated integral? Just use the Stream x that we previously defined: > ghci> erf = 0 :> exp (-x^2) > ghci> take 10 (toList erf) > [0.0,1.0,-0.0,-2.0,0.0,12.0,0.0,-120.0,0.0,1680.0] This is not even close to the tailor series of the erf function. This to me looks more like x' = exp(-x^2) which is entirely unrelated to the integral.

deleted 4y ago
[deleted]
texaslonghorn5 4y ago
The two are related because erf(x) is the integral of exp(-t^2) from 0 to x. So this is correct if you interpret coefficients in the basis {x^n/n!} and drop the scaling factor of 2/sqrt(pi). The Taylor series of exp has a 1 in the zeroth coefficient, so integrating from 0 to x we get a 1 in the 1st coefficient of erf. exp(-x^2) = 1 - x^2 + x^4/2 - x^6/6 + ... = 1/0! - 2x^2/2! + 12x^4/4! - 120x^6/6! + ... erf(x) = x/1! - 2x^3/3! + 12x^5/5! - 120x^7/7! + ...

snake_plissken 4y ago

Can someone better explain what the author means early on when they describe the Fundamental Theorem of Calculus as f(a + x) = f(a) + f'(t)dt|x -> a I thought that the integration yields the area under the curve from x to a? So how does the that relate to the output of f(a + x) which is not an area quantity?

cat_man 4y ago
I think it's a typo. \int_a^x f'(t)dt should be f(x) - f(a), and that would be the fundamental theorem of calculus (or some corollary of it). It looks like either f(a+x) should be f(x) or the upper limit x should be x+a. The rest of the post deals with Taylor series, which are often written with f(a+x) for the function being expanded, so that might have mixed the author up when writing that. They don't directly use that expression elsewhere as far as I can tell, so it doesn't affect the rest of the post. If it's f(x) = f(a) + \int_a^x f'(t)dt, then you can think of f(x) - f(a) as the area under the curve of f'(t) on the interval from a to x. For example, if f(x) = x, f'(x) = 1. The area of a function with height of 1 from a to x would be x - a, which is the same as f(x) - f(a) = x - a.
deleted 4y ago
[deleted]

gaze 4y ago

This seems a bit icky as you move away from the origin. There’s no promises for how fast the series converges and how well it does yea far out and so on and so forth. This reminds me a bit of harmonic balance which is a great, well behaved spectral technique