#set document(title: "1.7 Numerical methods: Euler’s method", author: "Jiří Lebl") #set page(width: 8.5in, height: auto, margin: 1in) #import "@preview/cetz:0.5.2" #set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en") #show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math")) #set par(justify: true, leading: 0.62em, spacing: 0.9em) #set enum(spacing: 1.1em) // room between list items so tall inline fractions don't collide #set list(spacing: 1.1em) #set table(stroke: 0.5pt + rgb("#c7ccd3")) #let BLUE = rgb("#183B6F") // brand navy — section bars + example/solution labels (white on navy 11.09:1) #let ORANGE = rgb("#A94509") // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text) #let RED = rgb("#DC2626") // brand error-600 #let GREEN = rgb("#059669") // brand success-600 (decoration only; small green text uses green-text #007942) #show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt, fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt, text(fill: white, weight: "bold", size: 19pt, it.body)) #show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt, fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt, text(fill: white, weight: "bold", size: 12pt, it.body)) #show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body) #show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body) #let examplebox(label, title, body) = block(width: 100%, breakable: true, fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt, above: 12pt, below: 12pt)[ #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt, text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]] #body] // rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade #let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint, stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt), above: 11pt, below: 11pt)[ #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)] #linebreak() #body] #let solutionbox(body) = block(above: 4pt, below: 8pt)[ #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body] #let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"), stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[ #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))] // Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a // consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the // print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a // degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero. #let _STD_W = 3.5 #let _WIDE_W = 5.6 #let _MAX_H = 3.4 #let _ASPECT_WIDE = 2.2 #let _UPSCALE_MAX = 1.15 #let stdfig(body) = context { let m = measure(body) let w = calc.max(m.width / 1in, 0.01) let h = calc.max(m.height / 1in, 0.01) let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W } let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX) align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body))) } #show figure: set block(breakable: false) #set figure(gap: 8pt) #show figure.caption: set text(size: 8.5pt, fill: rgb("#555")) == 1.7#h(0.6em)Numerical methods: Euler’s method At this point it may be good to first try the Lab II and/or Project II from the IODE website: www.math.uiuc.edu/iode/. As we said before, unless #math.equation(block: false, alt: "f open parenthesis x , y close parenthesis")[$f ( x , y )$] is of a special form, it is generally very hard if not impossible to get a nice formula for the solution of the problem #math.equation(block: true, alt: "y prime equals f open parenthesis x , y close parenthesis , y open parenthesis x sub 0 close parenthesis equals y sub 0")[$y^(′) = f ( x , y ) , y ( x_(0) ) = y_(0)$] If the equation can be solved in closed form, we should do that. But what if we have an equation that cannot be solved in closed form? What if we want to find the value of the solution at some particular #math.equation(block: false, alt: "x")[$x$]? Or perhaps we want to produce a graph of the solution to inspect the behavior. In this section we will learn about the basics of numerical approximation of solutions. The simplest method for approximating a solution is Euler's Method.#math.equation(block: false, alt: "to the power 1")[$1$]It works as follows: Take #math.equation(block: false, alt: "x sub 0")[$x_(0)$] and compute the slope #math.equation(block: false, alt: "k equals f open parenthesis x sub 0 , y sub 0 close parenthesis")[$k = f ( x_(0) , y_(0) )$]. The slope is the change in #math.equation(block: false, alt: "y")[$y$] per unit change in #math.equation(block: false, alt: "x")[$x$]. Follow the line for an interval of length #math.equation(block: false, alt: "h")[$h$] on the #math.equation(block: false, alt: "x")[$x$]-axis. Hence if #math.equation(block: false, alt: "y equals y sub 0")[$y = y_(0)$] at #math.equation(block: false, alt: "x sub 0")[$x_(0)$], then we say that #math.equation(block: false, alt: "y sub 1")[$y_(1)$] (the approximate value of #math.equation(block: false, alt: "y")[$y$] at #math.equation(block: false, alt: "x sub 1 equals x sub 0 plus h")[$x_(1) = x_(0) + h$]) is #math.equation(block: false, alt: "y sub 1 equals y sub 0 plus h k")[$y_(1) = y_(0) + h k$]. Rinse, repeat! Let #math.equation(block: false, alt: "k equals f open parenthesis x sub 1 , y sub 1 close parenthesis")[$k = f ( x_(1) , y_(1) )$], and then compute #math.equation(block: false, alt: "x sub 2 equals x sub 1 plus h")[$x_(2) = x_(1) + h$], and #math.equation(block: false, alt: "y sub 2 equals y sub 1 plus h k")[$y_(2) = y_(1) + h k$]. Now compute #math.equation(block: false, alt: "x sub 3")[$x_(3)$] and #math.equation(block: false, alt: "y sub 3")[$y_(3)$] using #math.equation(block: false, alt: "x sub 2")[$x_(2)$] and #math.equation(block: false, alt: "y sub 2")[$y_(2)$], etc. Consider the equation #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$], #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$], and #math.equation(block: false, alt: "h equals 1")[$h = 1$]. Then #math.equation(block: false, alt: "x sub 0 equals 0")[$x_(0) = 0$] and #math.equation(block: false, alt: "y sub 0 equals 1")[$y_(0) = 1$]. We compute #math.equation(block: true, alt: "x sub 1 equals x sub 0 plus h equals 0 plus 1 equals 1 , y sub 1 equals y sub 0 plus h f open parenthesis x sub 0 , y sub 0 close parenthesis equals 1 plus 1 times the fraction 1 over 3 equals the fraction 4 over 3 approximately equals 1.333 ,; x sub 2 equals x sub 1 plus h equals 1 plus 1 equals 2 , y sub 2 equals y sub 1 plus h f open parenthesis x sub 1 , y sub 1 close parenthesis equals the fraction 4 over 3 plus 1 times the fraction open parenthesis the fraction 4 over 3 close parenthesis squared over 3 equals the fraction 52 over 27 approximately equals 1.926 .")[$x_(1) = x_(0) + h = 0 + 1 = 1 , y_(1) = y_(0) + h f ( x_(0) , y_(0) ) = 1 + 1 · frac(1, 3) = frac(4, 3) ≈ 1.333 , \ x_(2) = x_(1) + h = 1 + 1 = 2 , y_(2) = y_(1) + h f ( x_(1) , y_(1) ) = frac(4, 3) + 1 · frac(attach(( frac(4, 3) ), t: 2), 3) = frac(52, 27) ≈ 1.926 .$] We then draw an approximate graph of the solution by connecting the points #math.equation(block: false, alt: "open parenthesis x sub 0 , y sub 0 close parenthesis")[$( x_(0) , y_(0) )$], #math.equation(block: false, alt: "open parenthesis x sub 1 , y sub 1 close parenthesis")[$( x_(1) , y_(1) )$], #math.equation(block: false, alt: "open parenthesis x sub 2 , y sub 2 close parenthesis")[$( x_(2) , y_(2) )$],…. For the first two steps of the method see Figure #math.equation(block: false, alt: "1")[$1$]. #figure(figph[Two copies of the same slope field, x from -1 to 3 and y from 0 to 3, slope depending only on height. A green Euler segment runs from (0, 1) to about (1, 1.33); in the right panel a steeper one continues to (2, 1.93).], alt: "Two copies of the same slope field, x from -1 to 3 and y from 0 to 3, slope depending only on height. A green Euler segment runs from (0, 1) to about (1, 1.33); in the right panel a steeper one continues to (2, 1.93).", caption: none) Figure #math.equation(block: false, alt: "1")[$1$]: First two steps of Euler's method with #math.equation(block: false, alt: "h equals 1")[$h = 1$] for the equation #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$] with initial conditions #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$]. More abstractly, for any #math.equation(block: false, alt: "i equals 0 , 1 , 2 , 3 , …")[$i = 0 , 1 , 2 , 3 , …$], we compute #math.equation(block: true, alt: "x sub i plus 1 equals x sub i plus h , y sub i plus 1 equals y sub i plus h f open parenthesis x sub i , y sub i close parenthesis .")[$x_(i + 1) = x_(i) + h , #h(2em) y_(i + 1) = y_(i) + h f ( x_(i) , y_(i) ) .$] The line segments we get are an approximate graph of the solution. Generally it is not exactly the solution. See Figure #math.equation(block: false, alt: "2")[$2$] for the plot of the real solution and the approximation. \# Euler's method on the book's own example: y' = y^2/3, y(0) = 1, find y(2). \# Pure Python - no imports, so it runs the instant the page loads. def euler(f, x0, y0, x\_end, h): x, y = x0, y0 while x \< x\_end - 1e-12: y += h \* f(x, y) x += h return y f = lambda x, y: y \* y / 3 \# the slope field this section plots exact = 3.0 \# y(2) = 3 for this initial-value problem print(f"{'h':\>9} {'Euler y(2)':\>12} {'error':\>10} {'ratio':\>7}") prev = None h = 1.0 for \_ in range(9): err = exact - euler(f, 0.0, 1.0, 2.0, h) ratio = f"{err / prev:7.3f}" if prev else " -" print(f"{h:9.5f} {exact - err:12.6f} {err:10.6f} {ratio}") prev, h = err, h / 2 print("\\nWatch the RATIO column, not the words: each row halves h, and the") print("error ratio creeps toward 0.5. That is what 'first order' means -") print("it only settles near 1/2 once h is small enough. Buying one more") print("decimal place costs about ten times the work.") #figure(figph[Slope field with x from -1 to 3 and y from 0 to 3 showing the exact solution as a smooth red curve and the two-step Euler approximation as a green polyline that stays below it, ending about one unit lower.], alt: "Slope field with x from -1 to 3 and y from 0 to 3 showing the exact solution as a smooth red curve and the two-step Euler approximation as a green polyline that stays below it, ending about one unit lower.", caption: [Figure #math.equation(block: false, alt: "2")[$2$]: Two steps of Euler’s method (step size 1) and the exact solution for the equation #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$] with initial conditions #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$].]) We continue with the equation #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$], #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$]. Let us try to approximate #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis")[$y ( 2 )$] using Euler’s method. In Figures #math.equation(block: false, alt: "1")[$1$] and #math.equation(block: false, alt: "2")[$2$] we have graphically approximated #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis")[$y ( 2 )$] with step size 1. With step size 1, we have #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis approximately equals 1.926")[$y ( 2 ) ≈ 1.926$]. The real answer is 3. We are approximately #math.equation(block: false, alt: "1.074")[$1.074$] off. Let us halve the step size. Computing #math.equation(block: false, alt: "y sub 4")[$y_(4)$] with #math.equation(block: false, alt: "h equals 0.5")[$h = 0.5$], we find that #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis approximately equals 2.209")[$y ( 2 ) ≈ 2.209$], so an error of about #math.equation(block: false, alt: "0.791")[$0.791$]. Table #math.equation(block: false, alt: "1")[$1$] gives the values computed for various parameters. #notebox("Your Turn", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[ Solve this equation exactly and show that #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis equals 3")[$y ( 2 ) = 3$]. ] The difference between the actual solution and the approximate solution we will call the error. We will usually talk about just the size of the error and we do not care much about its sign. The main point is, that we usually do not know the real solution, so we only have a vague understanding of the error. If we knew the error exactly …what is the point of doing the approximation? #figure(table( columns: 4, align: left, inset: 6pt, table.header([#math.equation(block: false, alt: "h")[$h$]], [Approximate #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis")[$y ( 2 )$]], [Error], [#math.equation(block: false, alt: "the fraction Error over Previous error")[$frac("Error", "Previous error")$]]), [1], [1.92593], [1.07407], [], [0.5], [2.20861], [0.79139], [0.73681], [0.25], [2.47250], [0.52751], [0.66656], [0.125], [2.68034], [0.31966], [0.60599], [0.0625], [2.82040], [0.17960], [0.56184], [0.03125], [2.90412], [0.09588], [0.53385], [0.015625], [2.95035], [0.04965], [0.51779], [0.0078125], [2.97472], [0.02528], [0.50913], )) Table #math.equation(block: false, alt: "1")[$1$]: Euler’s method approximation of #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis")[$y ( 2 )$] where of #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$], #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$]. We notice that except for the first few times, every time we halved the interval the error approximately halved. This halving of the error is a general feature of Euler’s method as it is a first order method. In the IODE Project II you are asked to implement a second order method. A second order method reduces the error to approximately one quarter every time we halve the interval (second order as #math.equation(block: false, alt: "the fraction 1 over 4 equals the fraction 1 over 2 x the fraction 1 over 2")[$frac(1, 4) = frac(1, 2) "x" frac(1, 2)$]). To get the error to be within 0.1 of the answer we had to already do 64 steps. To get it to within 0.01 we would have to halve another three or four times, meaning doing 512 to 1024 steps. That is quite a bit to do by hand. The improved Euler method from IODE Project II should quarter the error every time we halve the interval, so we would have to approximately do half as many “halvings” to get the same error. This reduction can be a big deal. With 10 halvings (starting at #math.equation(block: false, alt: "h equals 1")[$h = 1$]) we have 1024 steps, whereas with 5 halvings we only have to do 32 steps, assuming that the error was comparable to start with. A computer may not care about this difference for a problem this simple, but suppose each step would take a second to compute (the function may be substantially more difficult to compute than #math.equation(block: false, alt: "the fraction y squared over 3")[$frac(y^(2), 3)$]). Then the difference is 32 seconds versus about 17 minutes. Note: We are not being altogether fair, a second order method would probably double the time to do each step. Even so, it is 1 minute versus 17 minutes. Next, suppose that we have to repeat such a calculation for different parameters a thousand times. You get the idea. Note that in practice we do not know how large the error is! How do we know what is the right step size? Well, essentially we keep halving the interval, and if we are lucky, we can estimate the error from a few of these calculations and the assumption that the error goes down by a factor of one half each time (if we are using standard Euler). #notebox("Your Turn", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[ In the table above, suppose you do not know the error. Take the approximate values of the function in the last two lines, assume that the error goes down by a factor of 2. Can you estimate the error in the last time from this? Does it (approximately) agree with the table? Now do it for the first two rows. Does this agree with the table? ] Let us talk a little bit more about the example #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$], #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$]. Suppose that instead of the value #math.equation(block: false, alt: "y open parenthesis 2 close parenthesis")[$y ( 2 )$] we wish to find #math.equation(block: false, alt: "y open parenthesis 3 close parenthesis")[$y ( 3 )$]. The results of this effort are listed in Table #math.equation(block: false, alt: "2")[$2$] for successive halvings of #math.equation(block: false, alt: "h")[$h$]. What is going on here? Well, you should solve the equation exactly and you will notice that the solution does not exist at #math.equation(block: false, alt: "x equals 3")[$x = 3$]. In fact, the solution goes to infinity when you approach #math.equation(block: false, alt: "x equals 3")[$x = 3$]. #figure(table( columns: 2, align: left, inset: 6pt, table.header([#math.equation(block: false, alt: "h")[$h$]], [Approximate #math.equation(block: false, alt: "y open parenthesis 3 close parenthesis")[$y ( 3 )$]]), [1], [3.16232], [0.5], [4.54329], [0.25], [6.86079], [0.125], [10.80321], [0.0625], [17.59893], [0.03125], [29.46004], [0.015625], [50.40121], [0.0078125], [87.75769], )) Table #math.equation(block: false, alt: "2")[$2$]: Attempts to use Euler’s to approximate #math.equation(block: false, alt: "y open parenthesis 3 close parenthesis")[$y ( 3 )$] where of #math.equation(block: false, alt: "y prime equals the fraction y squared over 3")[$y^(′) = frac(y^(2), 3)$], #math.equation(block: false, alt: "y open parenthesis 0 close parenthesis equals 1")[$y ( 0 ) = 1$]. Another case when things can go bad is if the solution oscillates wildly near some point. Such an example is given in IODE Project II. The solution may exist at all points, but even a much better numerical method than Euler would need an insanely small step size to approximate the solution with reasonable precision. And computers might not be able to easily handle such a small step size. In real applications we would not use a simple method such as Euler’s. The simplest method that would probably be used in a real application is the standard Runge-Kutta method (see exercises). That is a fourth order method, meaning that if we halve the interval, the error generally goes down by a factor of 16 (it is fourth order as #math.equation(block: false, alt: "the fraction 1 over 16 equals the fraction 1 over 2 x the fraction 1 over 2 x the fraction 1 over 2 x the fraction 1 over 2")[$frac(1, 16) = frac(1, 2) "x" frac(1, 2) "x" frac(1, 2) "x" frac(1, 2)$]). Choosing the right method to use and the right step size can be very tricky. There are several competing factors to consider. \# Try it on an equation whose answer you know, so you can SEE the error. \# y' = y, y(0) = 1 -\> y(x) = e^x \# Edit f, y0 and x\_end and re-run. Nothing leaves your browser. def euler(f, x0, y0, x\_end, h): x, y = x0, y0 while x \< x\_end - 1e-12: y += h \* f(x, y) x += h return y f = lambda x, y: y \# \<- change me x0, y0 = 0.0, 1.0 x\_end = 1.0 e = 2.718281828459045 \# the true y(1) for this f print(f"{'steps':\>7} {'h':\>9} {'Euler':\>12} {'error':\>11}") for n in (1, 2, 4, 8, 16, 32, 64): h = (x\_end - x0) / n approx = euler(f, x0, y0, x\_end, h) print(f"{n:7d} {h:9.5f} {approx:12.7f} {e - approx:11.7f}") print("\\nThe error shrinks like h, not like h^2: doubling the steps") print("halves the error. That is what 'first order' costs you.") - Computational time: Each step takes computer time. Even if the function #math.equation(block: false, alt: "f")[$f$] is simple to compute, we do it many times over. Large step size means faster computation, but perhaps not the right precision. - Roundoff errors: Computers only compute with a certain number of significant digits. Errors introduced by rounding numbers off during our computations become noticeable when the step size becomes too small relative to the quantities we are working with. So reducing step size may in fact make errors worse. - Stability: Certain equations may be numerically unstable. What may happen is that the numbers never seem to stabilize no matter how many times we halve the interval. We may need a ridiculously small interval size, which may not be practical due to roundoff errors or computational time considerations. Such problems are sometimes called stiff. In the worst case, the numerical computations might be giving us bogus numbers that look like a correct answer. Just because the numbers have stabilized after successive halving, does not mean that we must have the right answer. We have seen just the beginnings of the challenges that appear in real applications. Numerical approximation of solutions to differential equations is an active research area for engineers and mathematicians. For example, the general purpose method used for the ODE solver in Matlab and Octave (as of this writing) is a method that appeared in the literature only in the 1980s. === Footnotes \[1\] Named after the Swiss mathematician Leonhard Paul Euler (1707–1783). The correct pronunciation of the name sounds more like “oiler.”