Euler’s Method with Python

Johanna M Debrecht

Page |1

Euler's Method with Python

Differential Equations

Lab Description

In this lab, we are going to explore Euler's method of solving first-order, initial value problem differential equations by writing a program in Python. You do not need to be an expert at Python, or even know the language yet to complete the lab. You also do not need to have python installed on your computer, but you will need access to the internet. All the code you need is contained in this document.

Remember, Euler's method uses tangent lines and the Linear Approximation function from calculus to approximate solutions to differential equations. The primary value of studying Euler's method is pedagogical, as it is a good introduction to the ideas used in numerical integration of differential equations (DEs). However, it has limited use because of the error that increases with each step.

Euler's Method Overview and Review

Given a first-order, initial value problem of the form:

=y f= ( x, y), y ( x0 ) y0 ,

we will build a basic program to approximate the solution using Euler's method. To begin, we need to decide over what interval x0 , xp we wish to find a solution. We will also need to decide how many points n within the interval that we wish to approximate. Thus, the lengths of our "steps" h will be given by

h = xp - x0 . n -1

Johanna M Debrecht

Page |2

This will generate n evenly spaced points x0 , x1,... xp , where xk = x0 + k h. Our goal is to

produce a table like the following, with the program filling in the y ( xk ) values (similar to Table

2.6.1 or Table 2.6.2 on p. 78 of the textbook).

xy x0 y0

x1 y ( x1 ) x2 y ( x2 )

( ) xp y xp

Recall that the accuracy of Euler's method improves if you make the "step" size h smaller; that is, if you take more points on the interval (increase n). The width of the interval will also affect the accuracy of the approximation. Larger intervals will be less accurate, so you would want to use more points on the interval.

We know that y ( x0 ) = y0 , which means that y( x0 ) = f ( x0 , y0 ). Using the linear approximation

function, (which is a tangent line approximation to the function), we can write

y ( x) y ( x0 ) + y( x0 )( x - x0 ) =y0 + f ( x0 , y0 )( x - x0 ), for all x close to x0.

Hence, we can approximate that

y ( x1 ) y1 = L ( x1 ) = y ( x0 ) + f ( x0 , y0 )( x1 - x0 ) = y0 + hf ( x0 , y0 ). Then, since y ( x1 ) y1, we can substitute into the DE to get

= y( x1 ) f ( x1, y ( x1 )) f ( x1, y1 ).

Therefore, the linear (tangent line) approximation function becomes

y ( x) y ( x1 ) + y( x1 )( x - x1 ) y1 + f ( x1, y1 )( x - x1 ), for all x close to x1.

Using the same process, we approximate that

y ( x2 ) y ( x1 ) + f ( x1, y1 )( x2 - x1 ) y2 =y1 + f ( x1, y1 )( x2 - x1 ) =y1 + hf ( x1, y1 ).

Johanna M Debrecht

Page |3

In this fashion, we can find approximations for y ( xk ) yk , where the yk are defined recursively

by the formula

yk+=1 yk + hf ( xk , yk ).

The table of values produced in this way generates a list of points that approximate the solution to the DE. This is Euler's method.

Coding Euler's Method Using Python: Part 1

Step 1 SageMath is a free open-source mathematics software system licensed under the GPL (General Public License). Through it, you can get free access to python, R (used in statistics), Octave, java, C++, fortran, SageMath, Julia, and others. Go to the following website and create an account (Collaborative Calculation in the Cloud).

Step 2 Click on the orange ball surrounded by the blue arcs on the upper left of your home page.

Click on "Create a New Project," give it a title, then click on "Create Project." When created, click on the title to open it.

Johanna M Debrecht

Step 3 At the top, in the banner, you will see the following. Click on "New."

Page |4

Click on "Sage worksheet" in the first row, upper left.

You should now see something like the following.

Step 4 We will be writing our code in the big empty space, starting where you see the horizontal blue line. We need to load some software packages, written by others, to assist us in our calculations

Johanna M Debrecht

Page |5

and graphing. For example, numpy is a library containing a vast number of routines for numerical calculation in python, and matplotlib contains packages for plotting and graphing. You will also need to put your name into the file. You can do that by enclosing it between triple quotation marks. Anything between triple quotation marks will be ignored by Python. (See Substep G for a note about quote marks.) For a single line comment, you can put a hashtag or # at the beginning of that one line.

Substep A Type in the following code:

Import Packages from External Libraries """Put Your Name Here""" import numpy as np from matplotlib import pyplot as plt

You should be able to copy/paste it in.

Substep B Now we need to enter some information about our initial value, the width of the interval, the number of data points we want approximated, and the size of our steps h. For the purposes of this lab, I am going to use the following DE as an example

y = 0.1 y + 0.4x2 , y (2) = 4.

Suppose we would like to approximate over the interval [2, 2.5] with n = 101. Then

( ) x0 = 2, y0 = 4, xp = 2.5, n = 101, and h = xp - x0 (n -1) = 0.005. Type in the following code to

enter this information.

Johanna M Debrecht

Page |6

Define Your Basic Values x0=2 y0=4 xp=2.5 n=101 step_h=(xp-x0)/(n-1)

Substep C In order to enter our x-values, we need to create a 101-tuple, or vector. The entries will be our x-

values, [2.0, 2.005, 2.01, 2.015,, 2.5]. Enter the following code to define this vector.

Create a Vector of the x-Values x=np.linspace(x0,xp,n)

Substep D Now that we have coded our x-values, we need to code how to generate the y-values. We will have to create a place to store the y-values. We will use an array (an ordered list with only one row or only one column) with n entries all set to 0.0. Enter the following code to create an array for storing the y-values.

Create an Array for the y-Values y=np.zeros([n])

Substep E We are now ready to code the part that uses Euler's method. We need to replace the entries in the y array (initially set to 0) with the approximations generated by Euler's method. Since this is an iterative process, we will use what is called a for loop in programming. This will cause the computer to repeat the steps inside the for loop a specified number of times.

The first line in the following box resets the first entry in the y array to the value of y0 = 4. (Remember that we set this in the code already with y0=4 in Substep B.)

Johanna M Debrecht

Page |7

The second line sets up a "dummy" variable to iterate for our for loop. Because the next line is indented (you must use a tab), then the variable will only be available for this level or one below it. When we move out of the tabbed area, it won't apply any more. The tab is critical; this is how python knows what belongs to what. Everything that is indented after the for will be considered part of the loop until the first line that is not indented.

The next line sets the kth entry of our y array to be the value of the (k - 1)th entry of the y array added to h times the approximated value of the DE's derivative. This line will input the approximated values of the solution to the DE into the y array, using Euler's method.

Create a for loop y[0]=y0 for k in range(1,n):

y[k]=y[k-1]+step_h*(0.1*(np.sqrt(y[k-1]))+0.4*(np.power(x[k-1],2)))

Substep F Since we have gone to all this trouble, we should definitely print out our results! We will use another for loop to print out our approximated values for the solution.

Print the Approximated Solution Values for k in range(n):

print(x[k],y[k])

Substep G Finally, we want to make a graph of our result. We can use the plot function from matplotlib to do this (in the second line we called it in as plt).

Caution! Depending on the font you are using, python may not like your single and/or double quote marks. If you see the error in red (below the code in Plot the Solution), that is what happened. You can type them directly into the worksheet. (You may not have noticed, but if they were incorrect, they would have shown up as red in the code.)

Johanna M Debrecht

Plot the Solution plt.plot(x,y,'o') plt.xlabel("Value of x") plt.ylabel("Value of y") plt.title("Approximate Solution with Euler's Method") plt.show()

Page |8

Let me guess--you fixed all the single and double quotes, but you still have an error. Did you notice that there is a single quote in "Euler's?" If that one is curved, python will give you an error message. You have two choices: replace it with the ascii straight single quote or put a u just before the text string in the double quotes. The u will cause it to ignore any non-ascii characters and to print them as is. (The Unicode for the ascii double quote is x22 then hit Alt + x. The Unicode for the ascii single quote is x27.) See the second option in the next box.

Option to Print Non-ascii Code plt.title(u"Approximate Solution with Euler's Method")

Step 5 It is finally time to run the code and see how precise your typing skills are!

Click the green arrow at the top left of the page that says "Run" to run your code.

 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download