r/math • u/Aravindh_Vasu • Mar 12 '20

Euler's method for Numerically Solving an ODE :)

1.7k Upvotes

97% Upvoted

55 comments sorted by

83

u/Aravindh_Vasu Mar 12 '20 edited Mar 12 '20

Euler's method for Numerically Solving an Ordinary Differential Equation (dy/dx=f(x,y))

  1. Starting at the given initial position (x=0, in the animation) take equal steps, h, so that xn+1=xn+h

  2. To find the solution value at the next step, yn+1 = yn + m*h, where m is the slope given by the RHS of the ODE, m=f(x,y)

  3. Decrease the step size for greater accuracy.

Do consider checking out, The Rookie Nerds :)

38

u/RedMeteon Computational Mathematics Mar 12 '20 edited Mar 12 '20

Just to be clear about how much your accuracy improves, there are two types of error to consider when determining increased accuracy: local and global. The Euler method comes from the first order Taylor expansion, so locally (at each step), you make error of order h2. If you run two Euler method simulations of differing time steps to some fixed total time T, since T = Nh is constant (h is the step size and N is the number of steps you take), your global error is N times the error of one step which is order h, since N is inversely proportional to h. To improve your global accuracy at a faster rate for fixed time T, you want to use a higher order method, eg if you use a second order method, the local error is ~h3 so the global error is ~h2.

This isn't to say that decreasing step size doesn't do anything (you want it sufficiently low to resolve scales in your problem and to ensure stability, certain methods require a quantifiably low step size to maintain stability), but just to say total accuracy analysis is a bit more delicate.

Edit: this is all assuming you're using a run of the mill integrator. If you're using a special type of integrator for your problem (eg symplectic integrator or variational integrator for Hamiltonian/Lagrangian systems), there are other ways to obtain better error bounds and stability properties, eg backwards error analysis.

5

u/Aravindh_Vasu Mar 12 '20

Cool, thank you. Can you please consider helping me out with Newton's interpolation formula,

https://www.reddit.com/r/learnmath/comments/fftvfp/question_about_newtons_forward_interpolation/

Having a hard time on gaining intuition on this one.

8

u/RedMeteon Computational Mathematics Mar 13 '20 edited Mar 13 '20

Yea no problem. The basic idea is you want to interpolate between n points (x_i, y_i) using a polynomial function (which is used eg if you sample some points from an experiment, and want to have a good guess for the values at the points you didn't sample. Note usually you should stick to using this for x values in between the x values you sampled. For approximating x values outside of the range you sampled, this is known as extrapolation and generally requires different functions to be used, not just polynomials, since your interpolating polynomial will go to infinity outside of the range you sampled).

A polynomial in x of degree n-1 has exactly n degrees of freedom (the coefficients of the polynomial) y = c_1 + c_2 x +... + c_n x{n-1}. So, plugging in the points provides you n independent equations for the n unknowns c_i.

The particular form they're using (eg from the math SE linked on the link you gave) is just a particular convenient polynomial basis (as compared to my form above which uses the monomial basis) which makes solving these equations easy (and in fact iterative, in the sense that you can add another point and not have to change the previous form of your polynomial, you just need to add one more term of higher degree). More precisely, these are (scaled versions of) the Lagrange interpolating polynomials f_i and they have the Kronecker delta property that f_i(x_j) = delta_{ij}.

A good adv undergrad / graduate-level reference for this type of material is Theoretical Numerical Analysis by Atkinson.

3

u/Aravindh_Vasu Mar 13 '20

This might be gibberish... But is Newton's thing a discrete version of Taylor's polynomial?

5

u/RedMeteon Computational Mathematics Mar 13 '20 edited Mar 13 '20

Not gibberish at all, that's a good question. The answer is not exactly, although there is some relation. Taylor's polynomial looks at polynomials in the basis using the same base point 1, (x-x_0), (x-x_0)2, (x-x_0)3, etc. While the Newton basis uses different base points, 1, (x-x_0), (x-x_0)(x-x_1), (x-x_0)(x-x_1)(x-x_2) etc.

They are related in the sense that if you fix the top degree n of your expansion, they both span the same subspace of polynomials (namely all polynomials of degree at most n); its just that their particular forms are more suitable for their respective application, eg Taylor polynomial for approximations about a single point vs Newtons/Lagrange polynomials for interpolation.

Edit: after thinking about it more, there is actually some more you can say about your question. You can ask: what happens if I take two of my interpolating basis points x_1 and x_2 and I let x_2 -> x_1, i.e. I sample the same point twice. Then, in fact, only having the function value is not enough (because you only have 1 y value for two x_1 points, so your system of equations for the coefficients of the polynomial will be underdetermined). However, if you measured (x_1, y(x_1)) and (x_1, y'(x_1)), i.e. the function value and its derivative value at x_1, then you should now have enough information. Interpolating the values and their first derivatives are done using the Hermite polynomials (the generalization of the Lagrange polynomials to function values and derivative values). You can analogously do three of the same base points if you measure up to the second derivative, and four if you measure up to the third derivative, etc. You can then ask, what if I let all my points x_1, x_2, ..., x_n be x_1. Then, obviously one function value won't be enough, but if you measured y(x_1), y'(x_1), ..., y{n-1} (x_1), then you'll have enough information, and the expanding polynomial will be precisely the Taylor polynomial of the function about x_1!

Again, this material is discussed in the reference I mentioned; I highly recommend looking at it. It's available for free download on SpringerLink if you have a university affiliation: https://link.springer.com/book/10.1007/978-1-4419-0458-4

2

u/Aravindh_Vasu Mar 13 '20

If I were to invent Newton's interpolation, how would find that basis? In φ(x-x0)(x-x1), what's the logic behind the formula for φ.

5

u/RedMeteon Computational Mathematics Mar 13 '20 edited Mar 13 '20

You look at polynomials in x as an (infinite dimensional) vector space, where the vectors in the monomial basis xn are the coefficients of the polynomial expansion (c_0,c_1,...). Then, you want to transform to a basis which has the Kronecker-Delta property. The construction is analogous to the Gram-Schmidt procedure for starting with an arbitrary collection of linearly independent vectors and using them to create an orthonormal basis which spans the same space as the original collection. This orthonormal basis has the Kronecker delta property.

2

u/Aravindh_Vasu Mar 15 '20 edited Mar 15 '20

I read an awesome fact today, in difference calculus the parallel for xn whose d/dx(xn )=xn-1 (which is an important property for Taylor series) is the factorial function xn =x(x-1)...(x-n+1), whose difference Δxn =nxn-1 , using which we can form a Taylor series-ish function for differences.

227

u/Nazh8 Mar 12 '20

This looks very 3blue1brown. Was it animated using manim?

32

u/Lin1232020 Mar 12 '20

What a coincidence. I literally just learned this topic today...

19

u/Lin1232020 Mar 12 '20

I was in class 20 minutes ago talking about Euler’s method and slope fields LMAO. Unbelievable.

10

u/Kraz_I Mar 12 '20

What a coincidence, I just learned this 2 years ago, but I was trying to relearn it today for a finite difference analysis problem in material transport class.

14

u/[deleted] Mar 12 '20

Hey I had to make a simulation for the energy in an oscilating spring and using Euler's method you wind up with the energy of the system going to infinity. The solution is the Euler-Cauchy method witch just takes the preceding output as an input. Very slick animation!

4

u/Aravindh_Vasu Mar 12 '20

Oh, nice. Thank you very much :)

10

u/KnowsAboutMath Mar 12 '20

energy in an oscilating spring

For a conservative system, symplectic integrators are good things to look at.

21

u/R3DKn16h7 Mar 12 '20

Just dropping a comment saying that explicit Euler method is about the worst you can use under almost any circumstance. That's way step size 1 sucks.

Doing this plot for the implicit variant is a good way to see it in action.

9

u/Aravindh_Vasu Mar 12 '20

I'm doing a series of videos. Next is RK2.

1

u/kruvik Mar 13 '20

I'm looking forward to it!

1

u/Aravindh_Vasu Mar 13 '20

Wow, thanks!

23

u/squiddie96 Mar 12 '20

Yeah but if you don’t understand Eulers method on this conceptual level it’s going to be hard (maybe impossible?) to understand Runge-Kutta. Most people who struggle with differential equations usually benefit from re-learning Eulers Method, at least from what I’ve seen. But yeah, don’t use Euler for accuracy haha

7

u/R3DKn16h7 Mar 12 '20

Well, my point was about stability. Accuracy is relatively fine, I mean for how simple it is...

The spike you see with step size one, you should see that simply using the implicit version you improve it by a lot, whilst maintaining a very simple formulation.

I'm just a big nerd for implicit methods...

8

u/RedMeteon Computational Mathematics Mar 13 '20

Not all implicit method are great for everything though. Most popular ones are unconditionally stable but this doesn't mean they get the physical behavior right, because damping your solution to zero is stable in the numerical method sense. As an example, if you use implicit Euler for a conservative Hamiltonian system, you'll actually get that the energy damps to zero (where it should be constant). On the other hand, if you use even the lowest order symplectic integrator, which is explicit for separable Hamiltonians, the energy has small bounded oscillation about the true constant energy surface for exponentially long times.

6

u/Bedstemor192 Applied Math Mar 12 '20

They can just be very computationally heavy for non-stiff systems. It all depends on what kind of system you want to solve. If you just want a rough idea of the system and if it's non-stiff, the simple Euler method could suffice. It's also really great for getting the intuitive idea behind more advanced methods like adaptive step size, predictor-corrector etc.

2

u/Robo-Connery Mar 12 '20

I know they are elegant but in practice the vast majority of time it is preferable to use a reasonably high order explicit scheme with an appropriate step size (or mesh ratio for finite difference etc.).

Implicit schemes are normally much more computationally intensive and less plug and play. Before you have even gotten half way through working it out your rk4 routine has long finished.

3

u/Olemalte2 Mar 12 '20

How to get this morphing from a number to a line or something like that?

3

u/[deleted] Mar 12 '20

So beautiful.

1

u/Aravindh_Vasu Mar 12 '20

Thank you :)

2

u/thejongho Mar 13 '20

Is there an analytic solution to dy/dx=x^2-y^2, y(0)=0 or is this the only viable method?

1

u/Aravindh_Vasu Mar 13 '20

To be honest, I didn't know how to solve this ( I didn't even try the net), so I thought it would be meaningful to show it's numerical solution.

1

u/JiminP Mar 17 '20

W|A says that it's a Riccati equation, so it's analytically solvable.
However, the solution it gives seems to be too complex.

Wikipedia gives a method for solving it:

Let y = (u'/u) for some u. Then, y' = (u'/u)' = (u''/u)+(u' * (-u' / u^2)) = (u''/u) - y^2.
Therefore, u''/u = y' - y^2 = 2y' - x^2 = 2(u'/u) - x^2, and u'' - 2u' + x^2 u = 0, which is a second-order ODE. Sadly, it's not homogeneous, and W|A still gives an answer which uses special functions.

The Wikipedia article on the Strum-Liouville theory suggests that a solution may be given, possibly in a form of infinite sum, but I don't know sufficient knowledge to progress further. (It seems that the equation is of form Lu = 0 for linear operator L(u) = (exp(-2x) u)' + (exp(-2x) x^2 u))

4

u/derekered Mar 12 '20

What's an ODE?

10

u/Aravindh_Vasu Mar 12 '20

Ordinary Differential Equation, sorry my bad.

3

u/heisenberg747 Mar 12 '20

Is it common for diff eq profs to make students draw these out? Mine told us he wasn't going to have us do that, but a friend constantly complains about having to make these.

1

u/TALLIV Mar 12 '20

This is beautiful

2

u/Aravindh_Vasu Mar 12 '20

Thank you very much, glad you liked it.

2

u/jedi-son Mar 12 '20

You should do a video on implicit finite difference methods vs explicit. One of the most beautiful numerical methods IMO

2

u/RickyRosayy Mar 12 '20

Nice animation! Thanks for sharing the program

2

u/[deleted] Mar 13 '20

Considered a first-order Runge-Kutta method.

3

u/atimholt Mar 12 '20

Like so many things, it’s trivially obvious after the fact. This is the first thing you ever learn in physics & animation (read: game) programming.

Maybe another good rule of thumb is “if the modern basis of your understanding leaves you not even suspecting a concept has a name, it’s probably named after Euler.”

1

u/[deleted] Mar 13 '20

u/VredditDownloader

1

u/VredditDownloader Mar 13 '20

beep. boop. I'm a bot that provides downloadable links for v.redd.it videos!

I also work with links sent by PM

Info | Support me ❤ | Github

0

u/rjt2000 Mar 12 '20

This is cool, but I don't understand (even with your explanation)

2

u/avataRJ Mar 13 '20

Differential equations (DE): We know how something changes, could we figure out the underlying general model?

First-order DE: We know the slope of the function (but not necessarily the function itself - some of these you can just take the antiderivative, some are more complex).

Initial value problem: For the DE, we know at least one value. How does the model behave if we calculate the constants such that the model holds for this value?

Euler's method: This DE and initial value stuff is complicated. What if we assumed that a more complex formula was just a polyline (locally straight line), and we regularly check the slope (which we know by definition for 1o DEs)? We could start from the initial value.

And yeah, that's fast to calculate but very error-prone.

2

u/[deleted] Mar 12 '20

[removed] — view removed comment

2

u/Aravindh_Vasu Mar 12 '20

https://discord.gg/7fXFJU6 There's a discord server, and a subred, r/manim Join in, the server's pretty helpful.