Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
3859 views
Tweet
Share
Share
sage.plot.graphics.Graphics.SHOW_OPTIONS['figsize']=5

#Euler's Method

Euler's Method is an algorithm used to construct approximate solutions to a differential equation of the form dydx=f(x,y)\displaystyle\frac{dy}{dx}=f(x,y) starting at an initial point (x0,y0)(x_0,y_0).

Since the differential equation dydx=f(x,y)\displaystyle\frac{dy}{dx}=f(x,y) tells us the slope of the tangent line at any point on the xy-plane, we can find the slope at (x0,y0)(x_0,y_0) and move along the tangent line some distance to a point (x1,y1)(x_1,y_1). Since the solution curve is close to its tangent line (as long as we're not too far from the point of tangency), the point (x1,y1)(x_1,y_1) is almost on the solution curve.

Now we repeat the process. Find the tangent line at (x1,y1)(x_1,y_1) using the differential equation, follow it for a short distance, and find a new point (x2,y2)(x_2,y_2). This point is also close to the solution curve.

Repeat the process as many times as you like.


The process is the same each time, so we can develop an iterated formula and automate the process.

Let's determine how to get from (xn,yn)(x_n,y_n) to (xn+1,yn+1)(x_{n+1},y_{n+1}).

First, we need to find the tangent line at (xn,yn)(x_n,y_n).

In general, the tangent line to a function y(x)y(x) at a the point x=ax=a has equation TL(x)=y(a)+y’(a)(x−a)TL(x)=y(a)+y’(a)(x-a).

In this case, the derivative is given by the differential equation, a=xna=x_n, and y(a)=yny(a)=y_n, so we have TL(x)=yn+f(xn,yn)(x−xn)TL(x)=y_n+f(x_n,y_n)(x-x_n).

Therefore, yn+1=TL(xn+1)=yn+f(xn,yn)(xn+1−xn)y_{n+1}=TL(x_{n+1})=y_n+f(x_n,y_n)(x_{n+1}-x_n)

We will move along the tangent line the same horizontal distance each step of the process. In other words, xn+1−xnx_{n+1}-x_n is a constant, called the “step size.” We will call this hh.

####Euler’s Formula

xn+1=xn+hx_{n+1}=x_n+hyn+1=yn+f(xn,yn)⋅hy_{n+1}=y_n+f(x_n,y_n)\cdot h

###Example 1

Use Euler's Method to approximate the solution curve to the differential equation dydx=x⋅y\displaystyle\frac{dy}{dx}=x\cdot y that passes through the point (0,1)(0,1). Plot the approximation for 0≤x≤20\le x \le 2.

We'll start with a small example by hand, and then we'll let the computer do the work.

We will use just 5 steps. That means the step size is h=2−05=25h=\frac{2-0}{5} = \frac{2}{5}.

We'll start with x0=0x_0=0 and y0=1y_0=1, and then we will calculate new x-coordinates with xn+1=xn+hx_{n+1}=x_n+h and new y-coordinates with yn+1=yn+xn⋅yn⋅hy_{n+1}=y_n+x_n\cdot y_n\cdot h.

So x1=0+25=25x_1=0+\frac{2}{5}=\frac{2}{5} and y1=1+0⋅1⋅25=1y_1=1+0\cdot 1\cdot \frac{2}{5}=1,

and x2=25+25=45x_2=\frac{2}{5}+\frac{2}{5}=\frac{4}{5} and y2=1+25⋅1⋅25=2925y_2=1+\frac{2}{5}\cdot 1\cdot \frac{2}{5}=\frac{29}{25},

and x3=45+25=65x_3=\frac{4}{5}+\frac{2}{5}=\frac{6}{5} and y3=2925+45⋅2925⋅25=957625y_3=\frac{29}{25}+\frac{4}{5}\cdot \frac{29}{25}\cdot \frac{2}{5}=\frac{957}{625},

and x4=65+25=85x_4=\frac{6}{5}+\frac{2}{5}=\frac{8}{5} and y4=957625+65⋅957625⋅25=3540915625y_4=\frac{957}{625}+\frac{6}{5}\cdot \frac{957}{625}\cdot \frac{2}{5}=\frac{35409}{15625},

and x5=85+25=2x_5=\frac{8}{5}+\frac{2}{5}=2 and y5=3540915625+85⋅3540915625⋅25=1451769390625y_5=\frac{35409}{15625}+\frac{8}{5}\cdot \frac{35409}{15625}\cdot \frac{2}{5}=\frac{1451769}{390625}.

Now let's plot these six points.

point([(0,1),(2/5,1),(4/5,29/25),(6/5,957/625),(8/5,35409/15625),(2,1451769/390625)])

The six points above are approximately on the solution curve. If we connect the points with straight lines, we will have an approximate solution curve.

Of course, just 5 steps is not enough to get a good approximation, so we'll use the computer with many more steps.

x,y=var('x y')
f(x,y)=x*y            #this is the function given by the differential equation
x0=0                  #initial value of x given in the problem
y0=1                  #initial value of y given in the problem
x_end=2               #the x-value you want to stop at
n=50                  #number of points to calculate
h=(x_end-x0)/n        #this calculates the step size for you
xlist=[x0];ylist=[y0] #we will use lists to keep track of all the x's and y's
for i in range(n):    #here we apply Euler's Formula
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]  #Note: RR converts to a floating-point number to avoid Sage taking too much time with exact values.
p=point([(xlist[i],ylist[i]) for i in range(n+1)])+line([[xlist[i],ylist[i]] for i in range(n+1)]);p

Here is a plot of our approximation (blue) along with the actual solution (red).

x=var('x')
y=function('y',x)
F(x)=desolve(derivative(y,x)==f(x,y),y,[x0,y0])
print 'The exact solution is y=',F(x)
p+plot(F(x),xmin=x0,xmax=x_end,color='red')
The exact solution is y= e^(1/2*x^2) 

︠b61faada-27b9-4103-bad8-bce14d7c91fbi︠
%md
We can make the approximation better by increasing $n$ (this decreases the step size).

If we want to plot the approximation past $x=2$, then we can change x_end. The approximation in this problem gets worse when we are farther away from our starting point.

The interactive box below allows us to change $n$ and x_end. Experiment with different values.

We can make the approximation better by increasing nn (this decreases the step size).

If we want to plot the approximation past x=2x=2, then we can change x_end. The approximation in this problem gets worse when we are farther away from our starting point.

The interactive box below allows us to change nn and x_end. Experiment with different values.

@interact
def _(x_end=2,n=50):
    x,y=var('x y')
    f(x,y)=x*y
    x0=0
    y0=1
    h=x_end/n
    xlist=[0];ylist=[1]
    for i in range(n):
        xlist+=[xlist[i]+h]
        ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
    p=point([(xlist[i],ylist[i]) for i in range(n+1)])+line([[xlist[i],ylist[i]] for i in range(n+1)])
    show(p+plot(e^(1/2*x^2),xmin=0,xmax=x_end,ymax=ylist[n],color='red'))
Interact: please open in CoCalc

###Example 2

Consider the initial value problem dydx=sin⁡(x)y,y(0)=0.5\displaystyle\frac{dy}{dx}=\sin(x)y,\quad y(0)=0.5.

Use Euler's Method to approximate y(10)y(10).

Modifying the routine from Example 1:

x,y=var('x y')
f(x,y)=sin(x)*y            #this is the function given by the differential equation
x0=0                  #initial value of x given in the problem
y0=0.5                  #initial value of y given in the problem
x_end=10               #the x-value you want to stop at
n=50                  #number of points to calculate
h=(x_end-x0)/n        #this calculates the step size for you
xlist=[x0];ylist=[y0] #we will use lists to keep track of all the x's and y's
for i in range(n):    #here we apply Euler's Formula
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]  #Note: RR converts to a floating-point number to avoid Sage taking too much time with exact values.
p=point([(xlist[i],ylist[i]) for i in range(n+1)])+line([[xlist[i],ylist[i]] for i in range(n+1)]);p
The exact solution of this equation using the method of separation of variables is y(x)=12e1−cos⁡(x).y(x)=\frac{1}{2}e^{1-\cos(x)}. What qualatative differences does the numerical approximation have that our exact solution does not? 
show(p+plot(0.5*e^(1-cos(x)),x,0,10,color='red'))
The approximate solution is decaying in amplitutue and has a "delay" in it's movement. Increasing the number of iterations helps with this, but doesn't eliminate it. 
x,y=var('x y')
f(x,y)=sin(x)*y            #this is the function given by the differential equation
x0=0                  #initial value of x given in the problem
y0=0.5                  #initial value of y given in the problem
x_end=10               #the x-value you want to stop at
n=200                  #number of points to calculate
h=(x_end-x0)/n        #this calculates the step size for you
xlist=[x0];ylist=[y0] #we will use lists to keep track of all the x's and y's
for i in range(n):    #here we apply Euler's Formula
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]  #Note: RR converts to a floating-point number to avoid Sage taking too much time with exact values.
p=point([(xlist[i],ylist[i]) for i in range(n+1)])+line([[xlist[i],ylist[i]] for i in range(n+1)]);show(p+plot(0.5*e^(1-cos(x)),x,0,10,color='red'))

@interact
def _(x_end=30,n=100):
    x,y=var('x y')
    f(x,y)=sin(x)*y
    x0=0
    y0=0.5
    h=x_end/n
    xlist=[x0];ylist=[y0]
    for i in range(n):
        xlist+=[xlist[i]+h]
        ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
    p=point([(xlist[i],ylist[i]) for i in range(n+1)])+line([[xlist[i],ylist[i]] for i in range(n+1)])
    show(p+plot(0.5*e^(1-cos(x)),xmin=0,xmax=x_end,ymax=1.5*max(ylist),color='red'))
Interact: please open in CoCalc

A Discussion About Error

Let’s look at how the error of Euler’s method changes depending on how many iterations we perform on a certain problem. In particular, how the error changes depending on how many iterations we utilize. We consider the same initial value problem as above, dydx=sin⁡(x)y\frac{dy}{dx}=\sin(x)y with y(0)=0.5y(0)=0.5. How much error is there at x=10x=10? We know the exact solution is y(x)=0.5e1−cos⁡(x)y(x)=0.5e^{1-\cos(x)}, and thus the value at x=10x=10 is y(10)=0.5e1−cos⁡(10)≈3.1453y(10)=0.5e^{1-\cos(10)}\approx 3.1453, We use Euler’s method to approximate the solution numerically using N=10N=10, 2020, 4040, 8080, 160160, and 320320 iterations.

x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=10        #10 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
3.078084072 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=20        #20 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
2.081033837 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=40        #40 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
1.275088250 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=80        #80 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
0.7209832181 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=160        #160 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
0.3859768012 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=320        #320 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
0.2000688603 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=640        #640 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
0.1018892070 
x,y=var('x y')
f(x,y)=sin(x)*y
x0=0
y0=0.5
x_end=10
n=1280        #1280 Iterations
h=(x_end-x0)/n
xlist=[x0];ylist=[y0]
for i in range(n):
    xlist+=[xlist[i]+h]
    ylist+=[ylist[i]+RR(f(xlist[i],ylist[i]))*h]
3.1453-ylist[n].n(digits=10)  #Actual Value minus approximate
0.05140414293

What do you notice about the size of the error as we keep doubling the number of iterations? (Double iterations -> Half the error) This is actually expected! Euler's method to approximate the solution of an initial value problem is called a first order method. Meaning that each doubling of iterations, halves the error. A second order method (and they are out there) would mean doubling the iterations quarters the error (1/2)2(1/2)^2. A common numerical method is called the Runge-Kutta method. This scheme is fourth order. Doubling the iterations decreases the error by a factor of 24=162^4=16. However, there is a tradeoff, these iterations are more invloved and take more time computationally.

###Interactive Euler's Method

The Sage code below will perform Euler's Method with interactive input. For your assignment, you may use this interact or copy and paste the code we used in the examples above.

x,y=var('x y')
@interact(auto_update=False)
def _(t1=text_control('<h3>Euler\'s Method</h3>',label=''),func=input_box(x*y,r'$\displaystyle\frac{dy}{dx}=$',width=25),x0=input_box(0,'Initial x value',width=10),y0=input_box(1,'Initial y value',width=10),x_end=input_box(1,'Final x value',width=10),n=input_box(100,'Number of steps',width=10),makeplot=checkbox(default=True,label='Plot a graph?'),makepred=checkbox(default=True,label='Approximate y value at final x?'),t2=text_control('If you uncheck both boxes, you will get no output.',label='')):
    f(x,y)=func
    h=(x_end-x0)/n
    a=[x0];b=[y0]
    for i in range(n):
        a+=[a[i]+h]
        b+=[b[i]+RR(f(a[i],b[i]))*h]
    if makeplot:
        p=point([(a[i],b[i]) for i in range(n+1)])+line([[a[i],b[i]] for i in range(n+1)])
        show(p)
        print ''
    if makepred:
        print 'The y value at x =',x_end,'is approximately',b[n].n(digits=10)
Interact: please open in CoCalc

#Euler's Method Assignment

###Question 1

Use Euler's Method to graph an approximate solution curve to dydx=x2−y2\displaystyle\frac{dy}{dx}=x^2-y^2 with initial value (1,2)(1,2). Graph on the interval from x=1x=1 to x=3x=3 and use n=20n=20 steps.

###Question 2

Consider the initial value problem dydx=6x2−3x2yy(0)=3\displaystyle\frac{dy}{dx}=6x^2-3x^2y\quad y(0)=3.

####Part a

Use Euler's Method with n=50n=50 steps to graph an approximate solution curve on the interval from x=0x=0 to x=1x=1.

####Part b

The exact solution of this differential equation is y=2+e−x3y=2+e^{-x^3}. Add a graph of this curve to your graph in part a.

###Question 3

Consider the initial value problem dydx=y2cos⁡(x),y(0)=0.5\displaystyle\frac{dy}{dx}=y^2\cos(x),\quad y(0)=0.5.

####Part a

Approximate y(π6)y\left(\frac{\pi}{6}\right) using Euler's Method with n=20n=20 (round to 10 decimal places).

####Part b

Use separation of variables to solve the differential equation for yy.

####Part c

Plug x=0x=0 and y=0.5y=0.5 into your solution and solve for the constant CC.

####Part d

Find the exact value of y(π6)y\left(\frac{\pi}{6}\right).

####Part e

Subtract your approximation (part a) from the exact value. This is the error.

####Part f

Approximate y(π6)y\left(\frac{\pi}{6}\right) using Euler's Method with n=50n=50 (round to 10 decimal places).

####Part g

Calculate the new error. [This should be smaller than the error in part e.]