Calculus 2, Lab 4: Differential Equations

Calculus 2, Lab 4: Differential Equations#

Today we will be using the Numpy and SciPy Python libraries in order to numerical approximate solutions to Ordinary Differential Equations using the Euler method. The session will build on the Jupyter notebooks which you have been using in class this week.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Slope fields#

The following code was introduced in class to demonstrate how to use Python to generate a slope field for a differential equation of the form \(y'=f(x,y)\).

def dfield(f):
    for j in x:
        for k in y:
            slope = f(k,j)
            domain = np.linspace(j-0.05,j+0.05,2)
            def fun(y1,x1):
                z = slope*(domain-x1)+y1
                return z
            plt.plot(domain,fun(k,j),solid_capstyle='projecting',solid_joinstyle='bevel')

    plt.title("Slope field for y'=f(x,y)")
    plt.grid(True)
    plt.show()

Problem 1#

Produce slope fields for the following differential equations:

(a) \(y' = y\sin(x)\)

(b) \(y' = x^2(1-y)\)

For each equation, you will need to:

  • define the function f(y,x) (remember that we have to write y,x and not x,y!)

  • define x and y using the np.linspace command

  • run the code above by putting dfield(f) as the last line in your code cell.

Exact solutions#

Recall that for a separable differential equation of the form \(y'=f(x)g(y)\), we can solve by rearranging this as \(\dfrac{1}{g(y)}\dfrac{dy}{dx}=f(x)\) and then integrating both sides with respect to \(x\).

Problem 2#

Solve the initial value problems below (the equations are the same as in problem 1).

(a) \(y'=y\sin(x)\), with \(y(0)=\pi/2\).

(b) \(y'=x^2(1-y)\), with \(y(0)=3\).

You can either solve them by hand, or use SymPy. The SymPy code is as follows:

import sympy as sy
sy.init_printing()
x = sy.symbols('x')
y = sy.Function('y')

yprime = sy.diff(y(x),x)
eq = sy.Eq(yprime,x-y(x))
eq

What this does:

  • defines x as a symbolic variable

  • defines y (symbolically) as a function of x

  • defines yprime as the derivative of y with respect to x

  • defines a differential equation “eq” of the form \(y'=x-y\). (You should replace x-y with an approprirate function.)

  • displays the differential equation

You can now ask the computer to solve the equation using the command sy.dsolve(eq). If the output looks messy, you can try wrapping this command in sy.simplify(). To specify that you want to solve for \(y\), you can use sy.dsolve(eq,y(x)), and to enter an intial condition \(y(0)=1\), try sy.dsolve(eq,y,ics={y(0):1}) as an optional argument (after the y(x) in sy.dsolve(eq,y)).

Euler’s Method#

  • Essentially a numerical counterpart to slope fields.

  • Recall the linear approximation $\(g(x_0+\Delta x) \approx g(x_0) + g'(x_0)\Delta x.\)$

  • Given \(y'=f(x,y)\), look for solution \(y=g(x)\) with initial condition \(g(x_0)=y_0\).

  • Notice \(g'(x_0) = f(x_0,y_0)\).

  • Proceed in “steps”, with fixed step size \(\Delta x\) as follows:

  1. \((x_0,y_0)\) is given.

  2. \(x_1 = x_0+\Delta x\) and \(y_1 = y_0+f(x_0,y_0)\Delta x\)

  3. \(x_2 = x_0+2\Delta x\) and \(y_2 = y_1+f(x_1,y_1)\Delta x\)

  4. \(x_3 = x_0+3\Delta x\) and \(y_3 = y_2+f(x_2,y_2)\Delta x\)

etc.

Code for Euler’s Method#

The following is borrowed from Patrick Walls at UBC. We define a function called odeEuler that implements the steps outlined above.

Input is an ODE of the form \(y'=f(y,x)\). Syntax is odeEuler(f,y0,x), where:

  • f is our function

  • y0 is the initial value for \(y\), corresponding to the value in x at index 0.

  • x is a NumPy array (the values to be used for the steps)

Before you can use it, you will need to define the function f(y,x), and the values for x (for example, x=np.linspace(0,2,100) will use 100 samples for x between 0 and 2).

def odeEuler(f,y0,x):
    y = np.zeros(len(x))
    y[0] = y0
    for n in range(0,len(x)-1):
        y[n+1] = y[n] + f(y[n],x[n])*(x[n+1] - x[n])
    return y

The output is a NumPy array containing the \(y\) values correpsonding to the \(x\) values in x.

If you want to plot these points, you can do the following. (You can copy this code into the empty cells below, and then edit as needed.)

def f(y,x):
    return 2*x-3*y #you need to replace 2*x-3*y with your function here
x = np.linspace(a,b,n) #replace a, b, and n with numbers
y0 = 0 #replace 0 with the initial value for y here
y = odeEuler(f,y0,x)
plt.plot(x,y)
plt.show()

Problem 3#

Apply Euler’s method to the two initial value problems from Problem 2. For each problem, you will need to define the function f(y,x) and the initial value y0 before running the odeEuler and plot commands. For each differential equation, use \(x\) values from 0 to 2 with 50 steps.

Checking the solution#

How well does our algorithm work? We can plot our values against the exact ones.

For each initial value problem done in problems 2 and 3: first define the function f(y,x) used by the odeEuler function, and then define a function g(x) that gives the exact solution. You can use your answers from Problem 2, but you will need to enter them as NumPy functions rather than SymPy functions.

# You will need to change the definitions of f(y,x), y0 and g(x) below to the ones for part (a)
# You may also need to define x as an np.linspace
def f(y,x):
    return 2*x-3*y
def g(x):
    return (2/3)*x-2/9+(2/9)*np.exp(-3*x)
y0 = 0
x = np.linspace(0,2,50)

y_exact = g(x)
y_euler = odeEuler(f,y0,x)
plt.plot(x,y_euler,'b.-',x,y_exact,'r-')
plt.legend(['Euler','Exact'])
plt.axis([0,2,0,2])
plt.grid(True)
plt.show()

If your graphs don’t fit well in the window, you can adjust the \(x\) and \(y\) values shown by changing the numbers in the plt.axis line.

Next, repeat the above for part (b).

# enter f(y,x) and g(x) for part (b) here

# copy,paste, and re-run the plotting code above here

Additional Challenge#

If you have time, see if you can produce a plot that shows the slope field, the exact solution, and the Euler’s method approximation, all in the same plot.

Contents