A lazy way to solve differential equations | Hacker News Reader

A lazy way to solve differential equations

iagoleal.com

204 pointscosmic_quanta28 comments

In Julia, this is implemented with the number type of TaylorSeries.jl https://github.com/JuliaDiff/TaylorSeries.jl and captured in the integration package 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_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_e... and 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 .

Hahaha, my first thought as I clicked the comments was "I wonder if Chris already has this in DifferentialEquations.jl" and here we are.

I knew I'd seen something like this in Julia. Fantastic to see it's already implemented!

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.

> 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/

[2]: https://www.chebfun.org/

[3]: https://www.chebfun.org/about/projects.html

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

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 (last paragraph of sec 7.1)

Following the breadcrumbs to the IVP treatment in chapter 10

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.

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.

What exactly do you mean by multivariate?

More than one variable :) x/y/z etc

There are two ways differential equations can involve more than one variable.

dx/dt = stuff

dy/dt = stuff

dt/dx = stuff

dt/dy = stuff

Require very different approaches.

I think they might be asking about PDEs instead of ODEs in this case when they mentioned x, y, z as being multivariable.

I wasn't actually gonna use them to solve a diff equation though to be honest :) I was just experimenting with denoising diffusion and instead of a u-net to infer the denoising step, I was thinking about using a multivariate polynomial. (Since if you express the output denoising as a linear combination of polynomials finding the denoiser just turns into linear algebra. Might be a big solve etc but inferene is way fast and it's nice to know that you have the global opimum.) I'm sitting at the code right this instant and was using x, y, xy, x^2y, xy^2, etc, but this thread inspired me to try out chebychev as a possible alternative. :)

Any orthogonal polynomial family can be combined into a simple tensor product basis. So in your case you could try T_1(x), T_1(y), T_1(x)T_1(y), T_2(x), T_2(x)T_1(y), and so on.

Depending on what kind of 2d: https://www.chebfun.org/docs/guide/guide12.html or https://www.chebfun.org/docs/guide/guide16.html

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

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

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

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.

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.

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! + ...

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.

This post is very possible without infinite streams. A for loop or a while loop would do the trick just as well.

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.

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?

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.

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

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

A top hit for more info on this topic, is this 2015 HN comment thread: https://news.ycombinator.com/item?id=9617083

Reply on news.ycombinator.com