# 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`

commandrun 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 variabledefines

`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:

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

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

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

\(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.