[HN Gopher] A lazy way to solve differential equations
       ___________________________________________________________________
        
       A lazy way to solve differential equations
        
       Author : cosmic_quanta
       Score  : 153 points
       Date   : 2022-07-22 12:52 UTC (10 hours ago)
        
 (HTM) web link (iagoleal.com)
 (TXT) w3m dump (iagoleal.com)
        
       | ChrisRackauckas wrote:
       | 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 .
        
         | helmholtz wrote:
         | Hahaha, my first thought as I clicked the comments was "I
         | wonder if Chris already has this in DifferentialEquations.jl"
         | and here we are.
        
       | freemint wrote:
       | > 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 wrote:
         | 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).
        
       | johnbcoughlin wrote:
       | 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 wrote:
         | > 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
        
         | ChrisRackauckas wrote:
         | > 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 wrote:
           | 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.
        
         | naillo wrote:
         | 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.
        
           | freemint wrote:
           | What exactly do you mean by multivariate?
        
             | naillo wrote:
             | More than one variable :) x/y/z etc
        
               | freemint wrote:
               | There are two ways differential equations can involve
               | more than one variable.
               | 
               | dx/dt = stuff
               | 
               | dy/dt = stuff
               | 
               | OR
               | 
               | dt/dx = stuff
               | 
               | dt/dy = stuff
               | 
               | Require very different approaches.
        
               | naillo wrote:
               | 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. :)
        
               | johnbcoughlin wrote:
               | 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.
        
               | wheelinsupial wrote:
               | I think they might be asking about PDEs instead of ODEs
               | in this case when they mentioned x, y, z as being
               | multivariable.
        
           | bigbillheck wrote:
           | Depending on what kind of 2d:
           | https://www.chebfun.org/docs/guide/guide12.html or
           | https://www.chebfun.org/docs/guide/guide16.html
        
       | snake_plissken wrote:
       | 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?
        
       | cosmic_quanta wrote:
       | 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 wrote:
         | This post is very possible without infinite streams. A for loop
         | or a while loop would do the trick just as well.
        
           | cosmic_quanta wrote:
           | 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 wrote:
        
       | lupire wrote:
       | 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
        
       | superb-owl wrote:
       | Nit: the FTC example has a typo - the integral should end at
       | `a+x`, not `a`
        
         | Blammar wrote:
         | 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 wrote:
       | 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.
        
         | texaslonghorn5 wrote:
         | 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! + ...
        
         | [deleted]
        
       ___________________________________________________________________
       (page generated 2022-07-22 23:00 UTC)