Direction Field 1st order ODE

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

matplotlib notebook

# ADJUST: This cell defines the right-hand side of the ODE
def RHS(T,Y):
    z = np.exp(T) - 2*Y
    return z

# Create the grid and calculate arrows at each point
# ADJUST: change the range of t and y here, inside linspace
T, Y = np.meshgrid(np.linspace(0,2,num=25), np.linspace(-1,2.5,num=35))
DT = np.ones(T.shape)
DY = RHS(T, Y)

# Scale the arrows to avoid too long and too short
# the choice of power (between 0 "don't scale" and 0.5 "all arrows equal length") to produce 
# the most aesthetically pleasing results depends on the differential equation
# ADJUST: pp
pp = 0.5
nn = (DT**2 + DY**2)**(pp)

# Plot the result
plt.figure();
# ADJUST: title and axes labels
plt.title('dy/dt = exp(t) - 2y'); plt.xlabel('t'); plt.ylabel('y');
plt.quiver(T, Y, DT/nn, DY/nn, angles='xy')
plt.grid()

## Also plot approximate solution calculated by advanced numerical methods (Runge-Kutta, chap 8)
sol = scipy.integrate.solve_ivp(RHS, [0,2], [1], max_step=0.02);
plt.plot(sol.t, sol.y[0,:], 'k-');

# For homework submission, you can save figure directly; uncomment the following
# plt.savefig('dirfield.pdf')

Euler Method

# ADJUST: This cell defines the right-hand side of the ODE, change as appropriate
def RHS(T,Y):
    z = (2-np.exp(T)) / (3+2*Y)
    return z

## Plot the direction field
# Create the grid and calculate arrows at each point
# ADJUST: change the range of t and y here, inside linspace
T, Y = np.meshgrid(np.linspace(0,1,num=26), np.linspace(-0,0.2,num=26))
DT = np.ones(T.shape)
DY = RHS(T, Y)
# Scale the arrows to avoid too long and too short. The choice of the exponent pp 
# (between 0 "don't scale" and 0.5 "all arrows equal length") 
# producing the most aesthetically pleasing results depends on the ODE.
# ADJUST: pp
pp = 0.5
nn = (DT**2 + DY**2)**(pp)

# Plot the result
# ADJUST: title and axes labels
plt.figure();
plt.title('dy/dt = (2-e^t) / (3+2y)'); plt.xlabel('t'); plt.ylabel('y');
plt.quiver(T, Y, DT/nn, DY/nn, angles='xy')
plt.grid()

## Euler method approximate solution ###########

# ADJUST: Euler step size
# h=0.1;
h = (T[0,1]-T[0,0])/0.2  # Step is dictated by the direction field grid step
# ADJUST: List of initial conditions for Y
Y0list = [ 0 ]

tt = np.arange(T.min(), T.max()+0.001, h)
for Y0 in Y0list:
    yy = np.zeros(tt.shape)
    yy[0] = Y0
    jmax = np.size(tt)
    for j in range(np.size(tt)-1):
        yy[j+1] = yy[j] + h*RHS(tt[j], yy[j])
        if yy[j+1] > Y.max() or yy[j+1] < Y.min():
            jmax = j+1
            break
    plt.plot(tt[:jmax],yy[:jmax], 'b-*')
    
## Also plot approximate solution calculated by advanced numerical methods (Runge-Kutta, chap 8)
#sol = scipy.integrate.solve_ivp(RHS, [0,1], [0], max_step=0.02);
#plt.plot(sol.t, sol.y[0,:], 'k-');

Direction Field and Implicit Solution

# ADJUST: This function defines the right-hand side of the ODE
def RHS(T,Y):
    z = (1 + 3*T**2) / (3*Y**2 - 6*Y)
    return z

# ADJUST: Initial condition
T0 = 0; Y0 = 1;

# ADJUST: This function defines our solution of the ODE in the *implicit* form F(t,y)=C
def Fsol(T,Y):
    z = Y**3 - 3*Y**2 - (T + T**3)
    return z

# ADJUST: The value of C for the solution passing through (T0, Y0)
C0 = -2;

# Create the grid and calculate arrows at each point
# ADJUST: change the range of t and y here, inside linspace
T, Y = np.meshgrid(np.linspace(-2,2,num=25), np.linspace(-1,2.5,num=35))
DT = np.ones(T.shape)
DY = RHS(T, Y)

# Scale the arrows to avoid too long and too short
# the choice of power (between 0 "don't scale" and 0.5 "all arrows equal length") to produce 
# the most aesthetically pleasing results depends on the differential equation
# ADJUST: pp
pp = 0.5
nn = (DT**2 + DY**2)**(pp)

# Plot the result
plt.figure();
plt.quiver(T, Y, DT/nn, DY/nn, angles='xy')
plt.contour(T,Y, Fsol(T,Y), levels=[C0], colors=['b','r'])
plt.plot(T0, Y0, '*', markersize=12)
plt.grid(); plt.xlabel('t'); plt.ylabel('y'); plt.title('ODE y\'=(1+3t^2)/(3y^2-6y)');

Euler vs Improved Euler method (chapter 8)

## Define the ODE
# ADJUST: This cell defines the right-hand side of the ODE, change as appropriate
def RHS(T,Y):
    z = (2-np.exp(T)) / (3+2*Y)
    return z

## Plot the direction field
# Create the grid and calculate arrows at each point
# ADJUST: change the range of t and y here, inside linspace
T, Y = np.meshgrid(np.linspace(0,1,num=25), np.linspace(-0,0.2,num=25))
DT = np.ones(T.shape)
DY = RHS(T, Y)
# Scale the arrows to avoid too long and too short. The choice of the exponent pp 
# (between 0 "don't scale" and 0.5 "all arrows equal length") 
# producing the most aesthetically pleasing results depends on the ODE.
# ADJUST: pp
pp = 0.5
nn = (DT**2 + DY**2)**(pp)

# Plot the result
plt.figure();
# ADJUST: title and axes labels
plt.title('dy/dt = (2-e^t)/(3+2y)'); plt.xlabel('t'); plt.ylabel('y');
plt.quiver(T, Y, DT/nn, DY/nn, angles='xy')
plt.grid()

## Euler method approximate solution ###########

# ADJUST: Euler step size
# h=0.1;
h = (T[0,1]-T[0,0])/0.2  # Step is dictated by the direction field grid step
# ADJUST: List of initial conditions for Y
Y0list = [ 0 ]

tt = np.arange(T.min(), T.max()+0.001, h)
for Y0 in Y0list:
    # Euler method
    yy = np.zeros(tt.shape)
    yy[0] = Y0
    jmax = np.size(tt)
    for j in range(np.size(tt)-1):
        yy[j+1] = yy[j] + h*RHS(tt[j], yy[j])
        if yy[j+1] > Y.max() or yy[j+1] < Y.min():
            jmax = j+1
            break
    plt.plot(tt[:jmax],yy[:jmax], 'b-*')
    
    # Improved Euler method
    yy = np.zeros(tt.shape)
    yy[0] = Y0
    jmax = np.size(tt)
    for j in range(np.size(tt)-1):
        k1 = RHS(tt[j], yy[j]);
        k2 = RHS(tt[j]+h, yy[j]+h*k1);
        yy[j+1] = yy[j] + h*(k1+k2)/2;
        if yy[j+1] > Y.max() or yy[j+1] < Y.min():
            jmax = j+1
            break
    plt.plot(tt[:jmax],yy[:jmax], 'r-*')
    
## Plot the high-precision numerical solution
sol = scipy.integrate.solve_ivp(RHS, [0,1], [0], max_step=0.02);
plt.plot(sol.t, sol.y[0,:], 'k-');