Lecture 1 Examples

By James Camacho

# Imports
import numpy as np
import scipy

%matplotlib inline
import matplotlib.pyplot as plt

# Approximate derivative of base^x:

print("For base 2:")
for dx in [1, 0.1, 0.01, 0.001]:
    print(f"dx = {dx} => (2^dx - 1) / dx = {(2**dx - 1) / dx}")
    
print("\nFor base 3:")
for dx in [1, 0.1, 0.01, 0.001]:
    print(f"dx = {dx} => (3^dx - 1) / dx = {(3**dx - 1) / dx}")

For base 2:
dx = 1 => (2^dx - 1) / dx = 1.0
dx = 0.1 => (2^dx - 1) / dx = 0.7177346253629313
dx = 0.01 => (2^dx - 1) / dx = 0.6955550056718884
dx = 0.001 => (2^dx - 1) / dx = 0.6933874625807412

For base 3:
dx = 1 => (3^dx - 1) / dx = 2.0
dx = 0.1 => (3^dx - 1) / dx = 1.1612317403390437
dx = 0.01 => (3^dx - 1) / dx = 1.104669193785357
dx = 0.001 => (3^dx - 1) / dx = 1.0992159842040383

# Calculating e:

print("Calculating e = lim_(n → ∞) (1 + 1/n)^n:")
for n in [10, 100, 1000, 10000]:
    print(f"For n = {n}, (1 + 1/n)^n = {(1+1/n)**n}")

print("\nCalculating e = 1/0! + 1/1! + 1/2! + ...")
e = 0
factorial = 1
for i in range(1, 100):
    e += 1 / factorial
    factorial *= i
print(f"e = {e:.4f}")

Calculating e = lim_(n → ∞) (1 + 1/n)^n:
For n = 10, (1 + 1/n)^n = 2.5937424601000023
For n = 100, (1 + 1/n)^n = 2.7048138294215285
For n = 1000, (1 + 1/n)^n = 2.7169239322355936
For n = 10000, (1 + 1/n)^n = 2.7181459268249255

Calculating e = 1/0! + 1/1! + 1/2! + ...
e = 2.7183

# Fixed point iteration for cos.
from math import cos

x = 1
while x != cos(x):
    x = cos(x)
print(f"Fixed point of cos is {x:.4f}.")

x = 1
def f(x):
    return (7 * x**3 - 3 * x**2 - 2) ** 0.2

while x != f(x):
    x = f(x)
print(f"Fixed point of (7x^3 - 3x^2 - 2)^(1/5) = {x:.3f}")

Fixed point of cos is 0.7391.
Fixed point of (7x^3 - 3x^2 - 2)^(1/5) = 2.362

# Newton's method for x^2 - 4

def f(x):
    return x**2 - 4

def f_prime(x):
    return 2*x

x = 3
while f(x) != 0:
    x -= f(x) / f_prime(x)
    print(x)
    
print(f"Found the solution x = {x}")

2.1666666666666665
2.0064102564102564
2.0000102400262145
2.000000000026214
2.0
Found the solution x = 2.0

# Euler's and trapezoidal methods

def y_true(t):
    return np.exp(t) - 1

def f(t, y):
    return y + 1

def euler(t, y, delta_t, t_end):
    points = [(t, y)]
    while t < t_end:
        y = y + delta_t * f(t, y)
        t += delta_t
        points.append((t, y))
    return points

def trapezoid(t, y, delta_t, t_end):
    points = [(t, y)]
    while t < t_end:
        y = ((2 + delta_t)*y + 2*delta_t) / (2 - delta_t)
        t += delta_t
        points.append((t, y))
    return points

t = np.linspace(0, 5)
plt.plot(t, y_true(t), label="True values", linewidth=5, color="yellow")

points = euler(0, 0, 0.05, 5)
t, y = zip(*points)
plt.plot(t, y, label="Euler's method")

points = trapezoid(0, 0, 0.05, 5)
t, y = zip(*points)
plt.plot(t, y, label="Trapezoidal method")

plt.title("y(t) = e^t - 1")
plt.xlabel("t")
plt.ylabel("y")
plt.legend()
plt.show()

# Euler local truncation error
deltas = np.linspace(0.0001, 0.1, 1000)
errors = []
for delta_t in deltas:
    t, y = 1, y_true(1)
    t, y = t+delta_t, y + delta_t*f(t, y)
    errors.append(abs(y - y_true(t)))

log_deltas = np.log(deltas)
log_errors = np.log(errors)
slope, _ = np.polyfit(log_deltas, log_errors, 1)

plt.title(f"LTE log-log plot, slope={slope:.2f}")
plt.plot(log_deltas, log_errors)
plt.show()

# Euler global runcation error
errors = []
for delta_t in deltas:
    t_end, y_end = euler(0, 0, delta_t, 5)[-1]
    errors.append(abs(y_end - y_true(t_end)))

log_deltas = np.log(deltas)
log_errors = np.log(errors)
slope, _ = np.polyfit(log_deltas, log_errors, 1)

plt.title(f"GTE log-log plot, slope={slope:.2f}")
plt.plot(log_deltas, log_errors)
plt.show()

# Adams Moulton coefficeints for s previous points.
s = 4

from scipy.special import factorial

n = np.arange(s+1)
k = np.arange(s+1)[:, None] # Turn [0, 1, ..., s] into a 1 x (s+1) matrix.

a_coeffs = np.r_[[0]*(s-1), -1, 1]
a_values = n**(k+1) / factorial(k+1) @ a_coeffs
B_equations = n**k / factorial(k)
b_coeffs = np.linalg.solve(B_equations, a_values)

print(f"a:\t{a_coeffs}\nb:\t{b_coeffs}")

def adams_bashforth(ts, ys, h, f):
    y = ys[-1]
    ys = np.r_[ys, [0]]
    ts = np.r_[ts, [ts[-1]+h]]
    for i in range(10): # Fixed point iteration
        a_piece = -a_coeffs[:-1] @ ys[:-1]
        hf = h * np.sum(b_coeffs @ f(ts, ys))
        ys[-1] = (a_piece + hf) / a_coeffs[-1]
    return ts[-1], ys[-1]

a:	[ 0  0  0 -1  1]
b:	[-0.02638889  0.14722222 -0.36666667  0.89722222  0.34861111]

# Example:
def y_true(t):
    return np.cos(t)**2

def f(t, y):
    return -2*y*np.tan(t)

h = 1e-3
ts = np.arange(s) * h
ys = y_true(ts)

t_end = 1

while ts[-1] < t_end:
    t, y = adams_bashforth(ts, ys, h, f)
    ts = np.roll(ts, -1)
    ys = np.roll(ys, -1)
    ts[-1] = t
    ys[-1] = y
    
print(f"Solved and true values at time {t:.3f}:\n{ys[-1]:.5f}, {y_true(t):.5f}")

Solved and true values at time 1.000:
0.29193, 0.29193

# Order of Adams Moulton method with s prior points is s+1
def error(f, y_true, t_start, dt, t_end):
    ts = t_start + dt * np.arange(s)
    ys = y_true(ts)
    
    while ts[-1] < t_end:
        t, y = adams_bashforth(ts, ys, dt, f)
        ts = np.roll(ts, -1)
        ys = np.roll(ys, -1)
        ts[-1] = t
        ys[-1] = y
    return abs(y_true(t) - y)

errors = []
deltas = np.linspace(0.001, 0.05)
for delta_x in deltas:
    e = error(f, y_true, 0, delta_x, 1)
    errors.append(e)

log_deltas = np.log(deltas)
log_errors = np.log(errors)
slope, intercept = np.polyfit(log_deltas, log_errors, 1)
delts = np.linspace(np.amin(log_deltas), np.amax(log_deltas))

plt.title(f"GTE log-log plot, slope={slope:.2f} (s={s})")
plt.plot(log_deltas, log_errors)
plt.plot(delts, slope * delts + intercept)
plt.show()