Pattern from Euler's Method on e^x

Art of Problem Solving

AoPS Online

College Math Topics in undergraduate and graduate studies

Topics in undergraduate and graduate studies

3 M G

BBookmark VNew Topic kLocked

College Math Topics in undergraduate and graduate studies

Topics in undergraduate and graduate studies

3 M G

BBookmark VNew Topic kLocked

w

2 users

TOPICS

No tags match your search

M

No tags match your search

M

Pattern from Euler's Method on e^x

G H J

Reveal topic tags

differential equations

inverse proportion

L

G H BBookmark kLocked kLocked NReply

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

1124 posts

fz0718

#1 h Jul 23, 2015, 4:45 PM • 2 Y

Y by Adventure10, Mango247

After reading the section in AoPS Calculus on Euler's Method, I was experimenting with a little python script to calculate $e^{0.3}$ with this method (as the book recommended), but I also tried varying the number of steps used.

import sys
import math
 
def diff_exp(x, y):
    '''A function representing a differential equation.
 
    This function represents y' = diff_equation(x, y), which
    is the differential equation for e^x.
    '''
    return y
 
def euler(start, val, target, equation, steps):
    increment = (target - start) / steps
    for step in range(steps):
        val = val + increment * equation(start, val)
        start += increment
    return val
 
def main():
    print("Steps|Absolute Error")
    print("-----|--------------")
    amt = 1
    for i in range(6):
        error = abs(euler(0, 1, 0.3, diff_exp, amt) - math.exp(0.3))
        print("10^%02d|%.12f" % (i, error))
        amt *= 10
 
if __name__ == '__main__':
    status = main()
    # sys.exit(status) # doesn't work on pywindow

import sys
import math

def diff_exp(x, y):
    '''A function representing a differential equation.

    This function represents y' = diff_equation(x, y), which
    is the differential equation for e^x.
    '''
    return y

def euler(start, val, target, equation, steps):
    increment = (target - start) / steps
    for step in range(steps):
        val = val + increment * equation(start, val)
        start += increment
    return val

def main():
    print("Steps|Absolute Error")
    print("-----|--------------")
    amt = 1
    for i in range(6):
        error = abs(euler(0, 1, 0.3, diff_exp, amt) - math.exp(0.3))
        print("10^%02d|%.12f" % (i, error))
        amt *= 10

if __name__ == '__main__':
    status = main()
    # sys.exit(status) # doesn't work on pywindow

If you run the code, you get the results from applying euler's method to calculate $e^{0.3}$ with $10^0$ , $10^1$ , $10^2$ , etc. steps.

But there appears to be a pattern in the errors. To see this pattern, here are a couple more rows of the output:

Steps|Absolute Error
-----|--------------
10^00|0.049858807576
10^01|0.005942428232
10^02|0.000606088209
10^03|0.000060730134
10^04|0.000006074229
10^05|0.000000607435
10^06|0.000000060744
10^07|0.000000006074

After $10^3$ steps, the error seems to decrease by almost exactly a factor of $10$ when the number of steps is increased by a factor of $10$ . Can anyone give an explanation for this seemingly inversely proportional behavior?

This post has been edited 1 time. Last edited by fz0718, Jul 23, 2015, 4:46 PM

Z K Y

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

196 posts

#2 h Jul 24, 2015, 3:31 AM • 4 Y

Y by fz0718, MSTang, Adventure10, Mango247

So the key in understanding why this occurs is to first understand what you're doing.

You're applying Euler's method, starting at $x=0$ , $y=1$ , and going to $x=.3$ in $n$ steps. At every step, you transform $f(x)$ into $f(x)+f(x)*.3/n$ , or in other words, each step takes $f(x)$ and replaces it with $f(x)*(1+.3/n)$ . Since we do this n times, the final value is $f(0)*(1+.3/n)^n=(1+.3/n)^n$ .

Now, when you find the error, you do the following computation: $exp(.3)-(1+.3/n)^n$ , and if we wanted to "normalize" this error, we could multiply by n and get $n(exp(.3)-(1+.3/n)^n)$ . At this point, the math becomes easier if we replace $n$ with $1/t$ :

$n(exp(.3)-(1+.3/n)^n=(exp(.3)-(1+.3t)^1/t))/t$

Our goal is to find what the limit is as n approaches infinity, or as t approaches 0. After applying L'Hopital's rule several times, including some logarithmic differentiation, you get a final answer: $.045exp(.3)=.0607436...$ . In general, when approximating the value of $e^a$ , the "normalized" error is $((a^2)*e^a)/2$ .

This post has been edited 1 time. Last edited by Grismor, Jul 24, 2015, 3:33 AM
Reason: Math typo

Z K Y

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

12096 posts

jmerry

#3 h Jul 24, 2015, 9:41 PM • 3 Y

Y by fz0718, Adventure10, Mango247

Actually, the specific differential equation has very little to do with this behavior, except that the exponential makes it easy to calculate exactly. If we use Euler's method to estimate the solution to a differential equation $y'=f(x,y)$ with steps of length $\epsilon$ , the truncation error at a single step is $\epsilon^2\cdot \frac{y''(x_0)}{2}$ , where $x_0$ is somewhere in that interval. To reach a fixed value of $x$ , we need a number of steps proportional to $\frac1{\epsilon}$ , for a total error proportional to $\epsilon$ .

But that's not all - we have errors compounding on errors, from the changes in $f$ when we input erroneous $y$ -values. How big are they? Well, first, let's assume a Lipschitz condition on $f$ , as in the uniqueness theorem. If $\left|f(x,y_1)-f(x,y_2)\right|\le L\left|y_1-y_2\right|$ , this will add at most $L\epsilon$ times the previous error to the error at each step. After $n$ steps, our first-step error is multiplied by at most $(1+L\epsilon)^n \le e^{L\cdot n\epsilon}=e^{L\delta x}$ . If we're targeting a fixed $x$ -value, that's a constant multiplier; compounding error makes our constant worse, but doesn't change the proportionality to step length.

If you don't have that Lipschitz condition, you can't expect Euler's method to work at all; that's the domain of such examples as non-unique solutions or solutions that blow up to $\infty$ at finite $x$ .

This post has been edited 1 time. Last edited by Kent Merryfield, Jul 25, 2015, 3:16 AM

Z K Y

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

Kent Merryfield

18574 posts

Kent Merryfield

#4 h Jul 25, 2015, 3:57 AM • 3 Y

Y by fz0718, Adventure10, Mango247

One amusing sideline is to try to solve the "differential equation" $y'=g(x).$ Of course, this is just an integration problem. Using Euler's method to approximate the solution is precisely the same thing as approximating the integral by left-endpoint Riemann sums. And the behavior of the error would be exactly the same as what you've documented.

jmerry has given the estimate that the error in a single step is proportional to $\epsilon^2$ where $\epsilon$ is the step size. Let me change letters of the alphabet and say the step size is $h,$ and the single-step error is $\frac12y''(\text{somewhere})\cdot h^2.$

In the conventional language of numerical analysis, we say that Euler's method is a first order method.

But wait a minute: isn't that an $h^2?$ Why are we calling it first order? It turns out that the proper definition of "local truncation error" is to take the (new) error produced by a single step of size $h$ , and then divide that error by $h$ . So in this case, we take the $O(h^2)$ that jmerry computed, divide it by $h,$ and get a local truncation error that is $O(h)$ -- in other words, first order. You can see the motivation for this definition in fz0718's post and code. He was aiming to cross a certain interval (in his case, to compute $y(0.3)$ starting from the initial condition $y(0)=1.$ ) To cross that interval with steps of size $h$ requires a number of steps that is proportional to $h^{-1},$ so we pile up that many new single-step errors. That justifies dividing by $h$ in the definition of local truncation error.

Any method for which the local truncation error tends to zero as $h\to 0$ is called consistent. A statement that the local truncation error tends to zero at a certain rate is a statement about the order of the method. And Euler's method is first order.

Any method that is not consistent is just formally wrong - you probably committed some sort of algebraic error when you designed it. But it turns out that that's not the only concern; you also have to worry about stability. Notice that above, I was careful to refer to the new error incurred at each step. But if you're partway across, your starting point for the step is itself in error. What happens as errors compound? I won't give a definition of stability - that's something to look up. It is possible to design a method that is consistent but unstable - but you'd have to use some kind of multi-step method to do that. A single-step method (like Euler) can't be unstable. Now what you really want of a method is for it to be convergent. That is, if we have $y(a)$ as a given initial condition and we want to compute $y(b)$ for some fixed $b,$ our estimate for $y(b)$ converges to the true $y(b)$ as $h\to 0.$ The big theorem: a method is convergent if and only if it is both consistent and stable.

One differential equations textbook that I know has a method that it calls "improved Euler." It turns out that this method also has another, more systematic name - but I'll leave that name out for now, so as not to spoil the fun. Here's the plan. Suppose we are trying to solve $y'=f(x,y),$ and we have a particular part-way-across estimated point $(x_n,y_n).$

Let $x_{n+1}=x_n+h.$
Form our first estimate of $y_{n+1}$ as in Euler. Call this $z$ and let $z=y_n+h\cdot f(x_n,y_n).$
Make a new change estimate based on this point and average the two to get the final estimate, which is to say:
Let $y_{n+1}=y_n+h\cdot(f(x_n,y_n)+f(x_{n+1},z))/2.$

(Of course, you should write the program so that you don't compute $f(x_n,y_n)$ twice.)

If you apply this method to the integration problem $y'=g(x),$ what well known method of integration do you get?

Now, I could whip up some numerical estimates myself easily enough, but I feel like having fz0718 practice his programming a little more. So here's the assignment (mostly intended for fz0718):

Program this method. Use the same problem as you did for Euler, and with the same goal - estimate $e^{0.3}.$ Use a variety of step sizes as before. (My one suggestion here is that cutting the step size by a factor of 10 each time may be a little too drastic - maybe you should use powers of 2 instead.) The errors should be going to zero, but how fast? From your numerical work, what does the order of this method appear to be?

This post has been edited 1 time. Last edited by Kent Merryfield, Jul 25, 2015, 3:58 AM

Z K Y

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

1124 posts

fz0718

#5 h Jul 25, 2015, 4:56 PM • 2 Y

Y by Adventure10, Mango247

Thank you all for your replies!

Code

import sys
import math
from decimal import Decimal
 
# thanks http://www-history.mcs.st-and.ac.uk/HistTopics/e_10000.html
E = Decimal('2.71828182845904523536028747135266249775'
            + '72470936999595749669676277240766303535')
 
def diff_exp(x, y):
    '''A function representing a differential equation.
 
    This function represents y' = diff_equation(x, y), which
    is the differential equation for e^x.
    '''
    return y
 
def euler(start, val, target, equation, steps):
    increment = (target - start) / steps
    for step in range(steps):
        val = val + increment * equation(start, val)
        start += increment
    return val
 
def euler_trapezoid(start, val, target, equation, steps):
    increment = (target - start) / steps
    for step in range(steps):
        result_1 = equation(start, val)
        result_2 = equation(start + increment, val + increment * result_1)
        val += increment * (result_1 + result_2) / 2
    return val
 
def main():
    print('Step|   Absolute Error   |Error Quotient')
    print('----|--------------------|--------------')
 
    prev = 0
    epow = E ** Decimal('0.3')
    for i in range(20):
        error = abs(euler_trapezoid(Decimal('0'), Decimal('1'), 
                    Decimal('0.3'), diff_exp, 1<<i) - epow)
        print('2^{:02d}|{:.18f}|{:.12f}'.format(i, error, prev / error))
        prev = error
 
if __name__ == '__main__':
    status = main()
    sys.exit(status)

Step|   Absolute Error   |Error Quotient
----|--------------------|--------------
2^00|0.004858807576003104|0.000000000000
2^01|0.001357245076003104|3.579904367980
2^02|0.000358884335331946|3.781845409184
2^03|0.000092279416044590|3.889104967445
2^04|0.000023396611524006|3.944135925406
2^05|0.000005890430914105|3.971969430620
2^06|0.000001477794475884|3.985960842479
2^07|0.000000370098646178|3.992974551910
2^08|0.000000092606020072|3.996485821228
2^09|0.000000023161681401|3.998242548524
2^10|0.000000005791692808|3.999121183939
2^11|0.000000001448082285|3.999560569414
2^12|0.000000000362040458|3.999780279071
2^13|0.000000000090512601|3.999890138127
2^14|0.000000000022628461|3.999945068711
2^15|0.000000000005657154|3.999972534268
2^16|0.000000000001414293|3.999986267112
2^17|0.000000000000353574|3.999993133549
2^18|0.000000000000088394|3.999996566786
2^19|0.000000000000022098|3.999998283501

The quotient between successive errors is 4 now, when the quotient between #steps is 2. It seems like the error is now proportional to the square of the number of steps, so it's a second-order method? Wow, it doesn't seem like going from left-endpoint riemann sums to trapezoid rule makes this much a difference!

I tried doing a similar thing with simpson's rule, but I still ended up with a degree 2 method. Do there exist, in general, arbitrarily high order methods?

Z K Y

The post below has been deleted. Click to close.

This post has been deleted. Click here to see post.

Kent Merryfield

18574 posts

Kent Merryfield

#6 h Jul 25, 2015, 7:42 PM • 3 Y

Y by fz0718, Adventure10, Mango247

Yes, "Improved Euler" is a second order method. The alternate name that I withheld earlier is "second order Runge-Kutta" or "RK2" for short. Applied to an integration problem, it would become the trapezoid rule, and yes, the trapezoid rule is second order. Going from left-endpoint to trapezoid would make precisely the same difference.

Simpson's rule is fourth order. The appropriate analogue for differential equations is fourth order Runge-Kutta, or RK4. RK4 is very much a workhorse of a method. Here's what it actually looks like (you wouldn't have guessed exactly this - you have to see it).

We're solving $y'=f(x,y)$ and we have reached point $(x_n,y_n).$ We're using step size $h.$
$\begin{align*}k_1 &:= f(x_n,y_n)\\ k_2 &:= f(x_n+h/2,y_n+h\cdot k_1/2)\\ k_3 &:= f(x_n+h/2,y_n+h\cdot k_2/2)\\ k_4 &:= f(x_n+h,y_n+h\cdot k_3)\\ x_{n+1}&:=x_n+h\\ y_{n+1}&:=y_n+h(k_1+2k_2+2k_3+k_4)/6\end{align*}$
One word of caution: in jmerry's post above, he hypothesized that $f(x,y)$ satisfied a Lipschitz condition in $y.$ That's strong enough to show that Euler is first order, but to show that some higher-order method is actually that higher order, you need for $f$ to have a certain number of bounded partial derivatives in $y$ as well.

One could construct methods of arbitrarily high order, but there are tradeoffs. The complexity of each step increases. If you were using a multi-step method, you have to worry about and analyze stability. And truncation error is always being played against roundoff error; the smaller the step size, the worse the problem with roundoff error. The main advantage of a higher order method is that it allows the use of larger step sizes, but beyond a certain point, you don't want to use a larger step size. RK4 hits a relative sweet spot among these tradeoffs.

Z K Y

N Quick Reply

G

H

=

© 2025 AoPS Incorporated

a