CoCalc Public Files14 - Euler's Method Assignment / Euler's Method Notes.sagewsOpen with one click!
Author: Aaron Tresham
Views : 47
License: Other -- explicitly state in your code
Compute Environment: Ubuntu 20.04 (Default)
This material was developed by Aaron Tresham at the University of Hawaii at Hilo and is Creative Commons License
licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.


  • Intro to Sage
  • Differential Equations

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 the point x=ax=a has equation TL(x)=y(a)+y(a)(xa)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)(xxn)TL(x)=y_n+f(x_n,y_n)(x-x_n).

Therefore, yn+1=TL(xn+1)=yn+f(xn,yn)(xn+1xn)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+1xnx_{n+1}-x_n is a constant, called the “step size.” We will call this hh.

In summary, we have the following:

Euler's Formula


yn+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=xy\displaystyle\frac{dy}{dx}=x\cdot y that passes through the point (0,1)(0,1). Plot the approximation for 0x20\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=205=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- and y-coordinates with the formulas xn+1=xn+hx_{n+1}=x_n+h yn+1=yn+xnynhy_{n+1}=y_n+x_n\cdot y_n\cdot h

x y
x0=0x_0=0 y0=1y_0=1
x1=0+25=25x_1=0+\frac{2}{5}=\frac{2}{5} y1=1+0125=1y_1=1+0\cdot 1\cdot \frac{2}{5}=1
x2=25+25=45x_2=\frac{2}{5}+\frac{2}{5}=\frac{4}{5} y2=1+25125=2925y_2=1+\frac{2}{5}\cdot 1\cdot \frac{2}{5}=\frac{29}{25}
x3=45+25=65x_3=\frac{4}{5}+\frac{2}{5}=\frac{6}{5} y3=2925+45292525=957625 y_3=\frac{29}{25}+\frac{4}{5}\cdot \frac{29}{25}\cdot \frac{2}{5}=\frac{957}{625}
x4=65+25=85x_4=\frac{6}{5}+\frac{2}{5}=\frac{8}{5} y4=957625+6595762525=3540915625 y_4=\frac{957}{625}+\frac{6}{5}\cdot \frac{957}{625}\cdot \frac{2}{5}=\frac{35409}{15625}
x5=85+25=2x_5=\frac{8}{5}+\frac{2}{5}=2 y5=3540915625+85354091562525=1451769390625 y_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.


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.

%var 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 p=point((x0,y0)) #we'll start keeping track of the points to graph 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=p+point((xlist[i+1],ylist[i+1]))+line([(xlist[i],ylist[i]),(xlist[i+1],ylist[i+1])]) show(p)

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


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. Of course, the approximation is going to get 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: please open in CoCalc

Example 2

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

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

I will copy and paste the formulas from above, skipping the plot:

%var y f(x,y)=y+x #dy/dx = y+x x0=0 #initial x-value = 0 y0=0 #initial y-value = 0 x_end=2 #we want to stop at x = 2 n=50 #we'll try 50 for now 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)] N(ylist[n]) #notice that ylist[n] is the last point calculated, in this case y(2)

We have found that y(2)4.1067y(2)\approx4.1067.

Let's try a higher value of n and see what happens.

%var y f(x,y)=y+x x0=0 y0=0 x_end=2 n=100 #we'll try n=100 this time 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)] N(ylist[n])

Now we have y(2)4.2446y(2)\approx4.2446.

Let's find the actual value. First, solve the differential equation.

y=function('y')(x) desolve(derivative(y,x)==y+x,y,[0,0])
-x + e^x - 1

Now plug in x=2x=2.

F(x)=-x + e^x - 1 #I'll call the solution F(x) F(2) N(F(2))
e^2 - 3 4.38905609893065

So y(2)=4.38905609893065y(2)=4.38905609893065.

Notice that increasing nn has gotten us closer to the actual answer. Let's increase nn one more time and see if we can get at least the first decimal place correct.

%var y f(x,y)=y+x x0=0 y0=0 x_end=2 n=500 #we'll try n=500 this time 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)] N(ylist[n])

Here is a summary of our results:

n Approximation
5050 4.106683346278314.10668334627831
100100 4.244646118252344.24464611825234
500500 4.359637175868974.35963717586897

The actual value is e234.38905609893065e^2-3\approx 4.38905609893065.