Copyright CRC Press 2023 - present CRC Press

Author: Dr Stephen Lynch National Teaching Fellow FIMA SFHEA

Stephen Lynch Homepage

Dear Readers: Feel free to contact me if you notice any errors. One of the advantages of showing solutions online is that they can be corrected.

Solutions for Section 1

Please note: Chapter 1 covers Python IDLE and you should run the commands below in IDLE.

I've used a jupyter notebook here so that the commands can be shown on the Web.

Chapter 1: The IDLE Integrated Development Learning Environment

# 1.1(a): Brackets.
2 * (3 - 4 * (5 - 8))

30

# 1.1(b): Exponential and trig function.
# Import all functions from the math library (module).
# Recall that pi/6 radians = 30 degrees.
from math import *
exp(sin(pi / 6))

1.648721270700128

# 1.1(c): Floor and ceiling functions.
floor(pi) - ceil(exp(1))

0

# 1.1(d): Floating point arithmetic.
3602879701896397 / 2**55

0.1

# 1.1(e): Fractions.
from fractions import Fraction
Fraction(4, 5) - Fraction(1, 7) * Fraction(2, 3)

Fraction(74, 105)

# 1.2(a): Generating data.
l1 = list(range(4 , 401 , 2))
print(l1)

[4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 136, 138, 140, 142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186, 188, 190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210, 212, 214, 216, 218, 220, 222, 224, 226, 228, 230, 232, 234, 236, 238, 240, 242, 244, 246, 248, 250, 252, 254, 256, 258, 260, 262, 264, 266, 268, 270, 272, 274, 276, 278, 280, 282, 284, 286, 288, 290, 292, 294, 296, 298, 300, 302, 304, 306, 308, 310, 312, 314, 316, 318, 320, 322, 324, 326, 328, 330, 332, 334, 336, 338, 340, 342, 344, 346, 348, 350, 352, 354, 356, 358, 360, 362, 364, 366, 368, 370, 372, 374, 376, 378, 380, 382, 384, 386, 388, 390, 392, 394, 396, 398, 400]

# 1.2(b): Generating data.
l2 = list(range(-9 , 200 , 4))
print(l2)

[-9, -5, -1, 3, 7, 11, 15, 19, 23, 27, 31, 35, 39, 43, 47, 51, 55, 59, 63, 67, 71, 75, 79, 83, 87, 91, 95, 99, 103, 107, 111, 115, 119, 123, 127, 131, 135, 139, 143, 147, 151, 155, 159, 163, 167, 171, 175, 179, 183, 187, 191, 195, 199]

# 1.2(c): Lists of lists.
A=[[1,2,3] , [4,5,6] , [7,8,9]]
A[1][2]

6

# 1.2(d): Compute the mean.
a = [10, 3, 4, 7, 8, 2, 5, 3, 4, 12]
mean_list = sum(a) / len(a)
print(mean_list)

5.8

# 1.3(a): Convert Fahrenheit to Centigrade.
def F2C():
    F = float(input("Enter temperature in degrees Fahrenheit: "))
    C = 5 * (F - 32) / 9
    print("Temperature in Centigrade is {:08.4f} C".format(C))
    
F2C()

Enter temperature in degrees Fahrenheit: 113
Temperature in Centigrade is 045.0000 C

# 1.3(b): Sum primes.
def sum_primes(n):
    sum_p = 0
    for n in range(2, n + 1):
        if all(n % i for i in range(2, n)):
            sum_p += n
    print('The sum of the primes up to {:,} is {:,}'.format(n, sum_p))

sum_primes(100)

The sum of the primes up to 100 is 1,060

# 1.3(c): Guess the number game.
import random # Import the random module.

num_guesses = 0
name = input('Hi! What is your name? ')
number = random.randint(1, 20) # A random integer between 1 and 20.
print('Welcome, {}! I am thinking of an integer between 1 and 20.'.format(name))

while num_guesses < 6:
    guess = int(input('Take a guess and type the integer? '))
    num_guesses += 1

    if guess < number:
        print('Your guess is too low.')
    if guess > number:
        print('Your guess is too high.')
    if guess == number:
        break

if guess == number:
    print('Well done {}! You guessed my number in {} guesses!'.format(name, num_guesses))
else:
    print('Sorry, you lose! The number I was thinking of was {}'.format(number))

Hi! What is your name? Ste
Welcome, Ste! I am thinking of an integer between 1 and 20.
Take a guess and type the integer? 4
Your guess is too low.
Take a guess and type the integer? 5
Your guess is too low.
Take a guess and type the integer? 15
Your guess is too high.
Take a guess and type the integer? 12
Your guess is too high.
Take a guess and type the integer? 10
Well done Ste! You guessed my number in 5 guesses!

# 1.3(d): Pythagorean triples.
import math
def pythagorean_triples(i):
    for b in range(i):
        for a in range(1, b):
            c = math.sqrt( a * a + b * b)
            n = 1
            if c - b == n:
                print(a, b, int(c))           
pythagorean_triples(101)

import math
def pythagorean_triples(i):
    for b in range(i):
        for a in range(1, b):
            c = math.sqrt( a * a + b * b)
            n = 3
            if c - b == n:
                print(a, b, int(c))
            
pythagorean_triples(201)

3 4 5
5 12 13
7 24 25
9 40 41
11 60 61
13 84 85
9 12 15
15 36 39
21 72 75
27 120 123
33 180 183

# 1.3(e): Words in a text.
import string
a_string = "One ring to rule them all, one ring to find them, one ring to bring them all, and in the darkness bind them."
# Ignore punctuation.
new_string = a_string.translate(str.maketrans('', '', string.punctuation))
print(new_string) 
# Split the words.
words = new_string.split(" ") 
long_words = []
count = 0
for j in range(len(words)):
  if len(words[j]) >= 5:
    count += 1
    long_words.append(words[j])
print("There are", count , "words with 5 letters or more.")
print("The words are:", long_words)

One ring to rule them all one ring to find them one ring to bring them all and in the darkness bind them
There are 2 words with 5 letters or more.
The words are: ['bring', 'darkness']

# 1.4: You must run this cell before the other turtle programs.
# These commands are not needed in IDLE.
!pip install ColabTurtlePlus
from ColabTurtlePlus.Turtle import *

Requirement already satisfied: ColabTurtlePlus in /usr/local/lib/python3.7/dist-packages (2.0.1)
Put clearscreen() as the first line in a cell (after the import command) to re-run turtle commands in the cell

# 1.4(a): Variation of the Cantor set.
initializeTurtle()
def cantor3(x , y , length):
  speed(13)
  if length >= 1:
    penup()
    pensize(2)
    pencolor("blue")
    setpos(x , y)
    pendown()
    fd(length)
    y -= 80
    cantor3(x , y , length / 5)
    cantor3(x + 2 * length / 5 , y , length / 5)
    cantor3(x + 4 * length / 5 , y , length / 5)
    penup()
    setpos(x , y + 80)
cantor3(-300 , 200 , 600)

# 1.4(b): The Koch square.
initializeTurtle()
pensize(1)
rt(90)
def KochSquare(length, level):
  speed(13)  # Fastest speed.
  for i in range(4):
    plot_side(length, level)
    rt(90)

def plot_side(length, level):
  if level==0:
    fd(length)
    return
  plot_side(length/3, level-1)
  lt(90)
  plot_side(length/3, level-1)
  rt(90)
  plot_side(length/3, level-1)
  rt(90)
  plot_side(length/3, level-1)
  lt(90)
  plot_side(length/3, level-1)
KochSquare(120 , 3)

# 1.4(c): A trifurcating tree.
initializeTurtle()
speed(13)
setheading(90)
penup()
setpos(0 , -200)
pendown()
def FractalTreeColor(length, level):
  pensize(length /10) # Thickness of lines.
  if length < 20:
    pencolor("green")
  else:
    pencolor("brown")
  if level > 0:
    fd(length)  # Forward
    rt(50)      # Right turn 50 degrees
    FractalTreeColor(length*0.7, level-1)
    lt(120)      # Left turn 120 degrees
    FractalTreeColor(length*0.5, level-1)
    rt(60)      # Right turn 60 degrees
    FractalTreeColor(length*0.5, level-1)
    rt(10)      # Right turn 10 degrees
    penup()
    bk(length)  # Backward
    pendown()
FractalTreeColor(200,6)

# 1.4(d): Sierpinski square.
initializeTurtle()
speed(13)
penup()
setpos(-200 , -200)
pendown()
def SierpinskiSquare(length, level):
  if level==0:
    return
  begin_fill() # Fill shape
  color("red")
  for i in range(4):
    SierpinskiSquare(length/3, level-1)
    fd(length/2)
    SierpinskiSquare(length/3, level-1)
    fd(length/1)
    lt(90) # Left turn 90 degrees
  end_fill()
SierpinskiSquare(200 , 3)

Chapter 2: Anaconda, Spyder and the Libraries NumPy, MatPlotLib and SymPy

# 2.1(a): An array. Vectorized computation.
# Sum the elements in each row.
import numpy as np
A = np.arange(100).reshape(10 , 10)
print("A = " , A)
A.sum(axis = 1)

A =  [[ 0  1  2  3  4  5  6  7  8  9]
 [10 11 12 13 14 15 16 17 18 19]
 [20 21 22 23 24 25 26 27 28 29]
 [30 31 32 33 34 35 36 37 38 39]
 [40 41 42 43 44 45 46 47 48 49]
 [50 51 52 53 54 55 56 57 58 59]
 [60 61 62 63 64 65 66 67 68 69]
 [70 71 72 73 74 75 76 77 78 79]
 [80 81 82 83 84 85 86 87 88 89]
 [90 91 92 93 94 95 96 97 98 99]]

array([ 45, 145, 245, 345, 445, 545, 645, 745, 845, 945])

# 2.1(b): Maxima of rows.
A.max(axis = 1)

array([ 9, 19, 29, 39, 49, 59, 69, 79, 89, 99])

# 2.1(c): Cumulative sum of columns.
A.cumsum(axis = 0)

array([[  0,   1,   2,   3,   4,   5,   6,   7,   8,   9],
       [ 10,  12,  14,  16,  18,  20,  22,  24,  26,  28],
       [ 30,  33,  36,  39,  42,  45,  48,  51,  54,  57],
       [ 60,  64,  68,  72,  76,  80,  84,  88,  92,  96],
       [100, 105, 110, 115, 120, 125, 130, 135, 140, 145],
       [150, 156, 162, 168, 174, 180, 186, 192, 198, 204],
       [210, 217, 224, 231, 238, 245, 252, 259, 266, 273],
       [280, 288, 296, 304, 312, 320, 328, 336, 344, 352],
       [360, 369, 378, 387, 396, 405, 414, 423, 432, 441],
       [450, 460, 470, 480, 490, 500, 510, 520, 530, 540]])

# 2.1(d): Finding an element in row 10 and column 5.
A[9][4]

94

# 2.1(e): Define a rank 3 tensor. A 2x2x2 array.
Tensor1 = np.array([[[1,2],[3,4]],[[5,6],[7,8]]])
print(Tensor1)
Tensor1.ndim

[[[1 2]
  [3 4]]

 [[5 6]
  [7 8]]]

3

# 2.2(a): Simple plot.
import matplotlib.pyplot as plt
x = np.linspace(-5,8,1000)
y = x**2 - 3 * x - 18
plt.plot(x , y)
plt.show()

# 2.2(b): Plot trigonometric function.
y = np.cos(2 * x)
plt.plot(x , y)
plt.xlabel("x")
plt.ylabel("y")
plt.show()

# 2.2(c): Plot.
y = np.sin(x)**2
plt.plot(x , y)
plt.xlabel("x")
plt.ylabel("y")
plt.title("$y=\sin^2(x)$")
plt.show()

# 2.2(d): Plot with point of inflection.
x = np.linspace(-2,2,1000)
y = 4 * x**3 - 3 * x**4
plt.ylim([-6, 2])
plt.plot(x , y)
plt.xlabel("x")
plt.ylabel("y")
plt.show()

# 2.2(e): Coshine graph.
x = np.linspace(-3 , 3 , 1000)
y = np.cosh(x)
plt.ylim([0, 10])
plt.plot(x , y)
plt.xlabel("x")
plt.ylabel("y")
plt.title("y=cosh(x)")
plt.show()

# 2.3(a): Factorize.
from sympy import *
x , y = symbols("x y")
factor(x**3 - y**3)

$\displaystyle \left(x - y\right) \left(x^{2} + x y + y^{2}\right)$

# 2.3(b): Solve.
# Solve x**2 - 7 * x - 30 = 0.
# With the solve function you do not have to include "=0".
solve(x**2 - 7 * x - 30 , x)

[-3, 10]

# 2.3(c): Partial fractions.
apart(3 * x / ((x-1) * (x+2) * (x-5)))

$\displaystyle - \frac{2}{7 \left(x + 2\right)} - \frac{1}{4 \left(x - 1\right)} + \frac{15}{28 \left(x - 5\right)}$

# 2.3(d): Simplify trig expressions.
trigsimp(sin(x)**4 - 2 * cos(x)**2 * sin(x)**2 + cos(x)**4)

$\displaystyle \frac{\cos{\left(4 x \right)}}{2} + \frac{1}{2}$

# 2.3(e): Expansion.
expand((y + x - 3) * (x**2 - y + 4)) 

$\displaystyle x^{3} + x^{2} y - 3 x^{2} - x y + 4 x - y^{2} + 7 y - 12$

# 2.4(a): Limits.
limit((1+1 / x)**x , x , oo)

$\displaystyle e$

# 2.4(b): Differentiation.
diff(3 * x**4 - 6 * x**3)

$\displaystyle 12 x^{3} - 18 x^{2}$

# 2.4(c): Higher order differentiation.
diff(3 * x**4 - 6 * x**3 , x , 3)

$\displaystyle 36 \left(2 x - 1\right)$

# 2.4(d): Indefinite integration.
print(integrate(x**2 - 2 * x - 3 , x), "+ c")

x**3/3 - x**2 - 3*x + c

# 2.4(e): Taylor series expansion.
(exp(x) * sin(x)).series(x , 1 , 10)

$\displaystyle e \sin{\left(1 \right)} + \left(x - 1\right) \left(e \cos{\left(1 \right)} + e \sin{\left(1 \right)}\right) + e \left(x - 1\right)^{2} \cos{\left(1 \right)} + \left(x - 1\right)^{3} \left(- \frac{e \sin{\left(1 \right)}}{3} + \frac{e \cos{\left(1 \right)}}{3}\right) - \frac{e \left(x - 1\right)^{4} \sin{\left(1 \right)}}{6} + \left(x - 1\right)^{5} \left(- \frac{e \sin{\left(1 \right)}}{30} - \frac{e \cos{\left(1 \right)}}{30}\right) - \frac{e \left(x - 1\right)^{6} \cos{\left(1 \right)}}{90} + \left(x - 1\right)^{7} \left(- \frac{e \cos{\left(1 \right)}}{630} + \frac{e \sin{\left(1 \right)}}{630}\right) + \frac{e \left(x - 1\right)^{8} \sin{\left(1 \right)}}{2520} + \left(x - 1\right)^{9} \left(\frac{e \cos{\left(1 \right)}}{22680} + \frac{e \sin{\left(1 \right)}}{22680}\right) + O\left(\left(x - 1\right)^{10}; x\rightarrow 1\right)$

# 2.5(a): Summation of an infinite series.
# Use oo for infinity.
n = symbols("n")
summation(1 / n , (n , 1 , oo) )

$\displaystyle \infty$

# 2.5(b): Solve linear simultaneous equations.
solve([x-y-2,x+y-1],[x,y])

{x: 3/2, y: -1/2}

# 2.5(c): Solve nonlinear simultaneous equations.
# You should also plot the curves to check the number of intersections.
solve([x**2-y-2,x+y-1],[x,y])

[(-1/2 + sqrt(13)/2, 3/2 - sqrt(13)/2), (-sqrt(13)/2 - 1/2, 3/2 + sqrt(13)/2)]

# 2.5(d): Arithmetic series.
summation(2 + 3 * (n-1) , (n , 1 , 20))

$\displaystyle 610$

# 2.5(e): Geometric series.
summation(2 * 3**n , (n , 1 , 20))

$\displaystyle 10460353200$

# 2.6(a): Matrix algebra.
A = Matrix([[-1,2,4],[0,3,2],[1,4,6]])
B = Matrix([[1,-1,1],[2,0,-1],[1,-1,1]])
3 * A - 5 * B

$\displaystyle \left[\begin{matrix}-8 & 11 & 7\\-10 & 9 & 11\\-2 & 17 & 13\end{matrix}\right]$

# 2.6(b): Matrix multiplication.
A * B

$\displaystyle \left[\begin{matrix}7 & -3 & 1\\8 & -2 & -1\\15 & -7 & 3\end{matrix}\right]$

# 2.6(c): Inverse matrix.
A.inv()

$\displaystyle \left[\begin{matrix}- \frac{5}{9} & - \frac{2}{9} & \frac{4}{9}\\- \frac{1}{9} & \frac{5}{9} & - \frac{1}{9}\\\frac{1}{6} & - \frac{1}{3} & \frac{1}{6}\end{matrix}\right]$

# 2.6(d): Access column 1.
(A**8).col(0)

$\displaystyle \left[\begin{matrix}771905\\575870\\1467293\end{matrix}\right]$

# 2.6(e): Eigenvalues and multiplicity.
B.eigenvals()

{1: 2, 0: 1}

# 2.6(e): You can also compute the eigenvectors.
# This was not requested in the question.
B.eigenvects()

[(0,
  1,
  [Matrix([
   [1/2],
   [3/2],
   [  1]])]),
 (1,
  2,
  [Matrix([
   [1],
   [1],
   [1]])])]

# 2.7(a): Complex number arithmetic.
# Generally, mathematicians use i for sqrt(-1) and engineers use j.
z1 , z2 = 1 - 4j, 5 + 6j 
3 * z1 - 5 * z2

(-22-42j)

# 2.7(b): Modulus of a complex number.
abs(z1 * z2)

32.202484376209235

# 2.7(c): Complex expand.
(z1 * exp(z2)).expand(complex=True)

$\displaystyle -23.3738416414842 - 611.476558862752 i$

# 2.7(d): Complex expand.
(sin(z1)).expand(complex=True)

$\displaystyle 22.9790855778861 - 14.7448051885587 i$

# 2.7(e): Convert to polar coordinates.
import cmath
cmath.polar(z2) 

(7.810249675906654, 0.8760580505981934)

Chapter 3: Jupyter Notebooks and Google Colab

# 3.1(a): Solve dy/dx=-y/x.
from sympy import *
x = symbols("x")
y = symbols("y", cls = Function)
ODE1 = Eq(y(x).diff(x) , - y(x) / x)
sol1 = dsolve(ODE1 , y(x))
print(sol1)

Eq(y(x), C1/x)

# 3.1(b): Solve dy/dx=y/x^2.
ODE2 = Eq(y(x).diff(x) , y(x) / x**2)
sol2 = dsolve(ODE2 , y(x))
print(sol2)

Eq(y(x), C1*exp(-1/x))

# 3.1(c): Solve dx/dt+x^2=1.
t = symbols("t")
x = symbols("x", cls = Function)
ODE3 = Eq(x(t).diff(t) , - x(t)**2 + 1)
sol3 = dsolve(ODE3 , x(t))
print(sol3)

Eq(x(t), -1/tanh(C1 - t))

# 3.1(d): Solve dx/dt+x=sin(t).
ODE4 = Eq(x(t).diff(t) , - x(t) + sin(t))
sol4 = dsolve(ODE4 , x(t))
print(sol4)

Eq(x(t), C1*exp(-t) + sin(t)/2 - cos(t)/2)

# 3.1(e): Solve y"+5y'+6y=10sin(t).
y = symbols("y", cls = Function)
ODE5 = Eq(y(t).diff(t,t) + 5 * y(t).diff(t) + 6 * y(t) , 10 * sin(t))
sol5 = dsolve(ODE5 , y(t))
print(sol5)

Eq(y(t), C1*exp(-3*t) + C2*exp(-2*t) + sin(t) - cos(t))

# 3.2: Animation.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

# Set up figure.
fig, ax = plt.subplots()
plt.title('Animation')
plt.close()
# Set domain and range.
ax.set_xlim(( -5, 5))
ax.set_ylim((- 5, 5))
# Set line width.
line, = ax.plot([], [], lw=2)
# Initialization function: plot the background of each frame.
def init():
    line.set_data([], [])
    return (line,)
# Animation function. This is called sequentially.  
def animate(n):
    x = np.linspace(-5 , 5 , 100)
    y = x**3 + (-4 + n * 0.1) * x + 1
    line.set_data(x, y)
    return (line,)
# Animate, interval sets the speed of animation.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                             frames=61, interval=100, blit=True)
# Note: below is the part which makes it work on Colab.
rc('animation', html='jshtml')
anim

Once Loop Reflect

# 3.3: Animation of parametric plot.
# Set up figure.
fig, ax = plt.subplots()
plt.title('Animation of a Parametric Curve')
plt.close()
# Set domain and range.
ax.set_xlim(( -1.2, 1.2))
ax.set_ylim((- 1.2, 1.2))
# Set line width.
line, = ax.plot([], [], lw=2)
# Initialization function: plot the background of each frame.
def init():
    line.set_data([], [])
    return (line,)
# Animation function. This is called sequentially.  
def animate(n):
    t = np.linspace(0, 2 * np.pi, 100)
    x = np.cos(t)
    y = np.sin(0.1 * n * t)
    line.set_data(x, y)
    return (line,)
# Animate, interval sets the speed of animation.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                             frames=41, interval=100, blit=True)
# Note: below is the part which makes it work on Colab.
rc('animation', html='jshtml')
anim

Once Loop Reflect

# 3.4: Interactive plot.
# Copy the code into your jupyter notebook.
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

from ipywidgets import interactive
import matplotlib.pyplot as plt
import numpy as np
def f(A, B, C):
    plt.figure(2)
    x = np.linspace(-4, 4, num=1000)
    plt.plot(x, A * np.exp(B * x) + C)
    plt.ylim(-5, 5)
    plt.show()
interactive_plot = interactive(f, A=(-2, 2 , 0.1), B = (-1, 1 , 0.1), \
                               C = (-3, 3 , 0.1))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

# 3.5: Interactive plot.
# Copy the code into your jupyter notebook.
from ipywidgets import interactive
import matplotlib.pyplot as plt
import numpy as np
def f(A , B , C , D):
    plt.figure(2)
    t = np.linspace(0, 8 * np.pi, num=1000)
    x = A * np.sin(B * t)
    y = C * np.cos(D * t)
    plt.plot(x, y)
    plt.ylim(-3, 3)
    plt.show()
interactive_plot = interactive(f, A=(-2, 2 , 0.1), B = (0, 2*np.pi , 0.1), \
                               C = (-2, 2 , 0.1) , D = (0, 2*np.pi , 0.1))
output = interactive_plot.children[-1]
output.layout.height = '350px'
interactive_plot

Chapter 4: Python for AS-Level Mathematics

# 4.1(a): Expansion.
expand((2 - sqrt(5))**4)

$\displaystyle 161 - 72 \sqrt{5}$

# 4.1(b): Solve.
# Solve 2 * x**2 - 3 * x - 2 = 0.
x , y = symbols("x y")
solve(2 * x**2 - 3 * x - 2 , x)

[-1/2, 2]

# 4.1(c): Solve.
solve([y - x - 1 , x**2 + y**2 - 2], [x , y])

[(-1/2 + sqrt(3)/2, 1/2 + sqrt(3)/2), (-sqrt(3)/2 - 1/2, 1/2 - sqrt(3)/2)]

# 4.1(d): Implicit plots.
fig = plt.figure(figsize = (6, 6))
x, y = np.mgrid[0:10:100j,0:10:100j]
z = x**2 + y**2 - 10 * x - 10 * y
# You can plot multiple contours if desired.
plt.contour(x, y, z, levels = [-45])
x = np.linspace(0, 10, 100)
y = -2 * x + 10
plt.plot(x, y)
plt.xlim(0 , 10)
plt.ylim(0 , 10)
plt.show()

# 4.1(e): Plots.
t = np.linspace(-2 , 2 , 1000)
y = 2 * np.sin(2 * np.pi * t)
plt.plot(t , y)
plt.show()

[<matplotlib.lines.Line2D at 0x7fc0ebe68f10>]

# 4.2(a): Factorize.
x = symbols("x")
factor(x**4 - 2 * x**3 - 6 * x**2 + 6 * x + 9)

$\displaystyle \left(x - 3\right) \left(x + 1\right) \left(x^{2} - 3\right)$

# 4.2(b): Two plots in one figure.
x = np.linspace(- 2 * np.pi , 2 * np.pi , 1000)
y1 , y2 = np.cos(x) , 1 - 2 * np.cos(x - np.pi)
plt.plot(x , y1 , x , y2)
plt.show()

[<matplotlib.lines.Line2D at 0x7fc0ec16f190>,
 <matplotlib.lines.Line2D at 0x7fc0ec192f90>]

# 4.2(c): Expansion.
x = symbols("x")
expand((1 - 2* x)**5)

$\displaystyle - 32 x^{5} + 80 x^{4} - 80 x^{3} + 40 x^{2} - 10 x + 1$

# 4.2(d): Differentiation.
print("y'=", diff(4 * x**8 - 3 * x**2 , x))
print('y"=', diff(4 * x**8 - 3 * x**2 , x , 2))

y'= 32*x**7 - 6*x
y"= 2*(112*x**6 - 3)

# 4.2(e): Integration.
integrate(3 * x**0.5 - x**1.5 , (x , 0 ,  3))

$\displaystyle 4.15692193816531$

# 4.3(a): Vectors in 2D.
u , v = np.array([2 , -3]) , np.array([1 , 4])
5 * u - 6 * v

array([  4, -39])

# 4.3(b): Solving power equations.
solve(5**(2 * x + 3) - 10 , x)

[(log(2/25)/2 + I*pi)/log(5), -1 + log(2)/(2*log(5))]

# 4.3(c): Generating data.
import math , random
a = list(range(1 , 500))
data1 = random.sample(a , 20)
print(data1)

[305, 491, 137, 464, 130, 186, 38, 69, 68, 390, 151, 278, 388, 371, 32, 320, 256, 60, 474, 197]

# 4.3(d): Basic statistics.
import statistics as stats
data1 = [2, 3, 5, 67, 46, 34, 34, 18, 4, 54, 24, 34, 56, 35, 66]
mean_data1 = stats.mean(data1)
print("mean_data1 =", mean_data1)
median_data1 = stats.median(data1)
print("median_data1 =", median_data1)
mode_data1 = stats.mode(data1)
print("mode_data1 =", mode_data1)
var_data1 = stats.variance(data1)
print("var_data1 =", var_data1)
stdev_data1 = stats.stdev(data1)
print("stdev_data1 =", stdev_data1)

mean_data1 = 32.13333333333333
median_data1 = 34
mode_data1 = 34
var_data1 = 512.2666666666667
stdev_data1 = 22.63330878741919

# 4.3(e): Venn diagrams.
plt.figure(figsize=(4,4))
from matplotlib_venn import venn3
fig, ax = plt.subplots(facecolor="lightgray")
plt.rcParams["font.size"] = "15"
ax.axis([0, 1, 0, 1])
v = venn3(subsets=(4,3,5,7,0,0,0),set_labels = ('Cycle', 'Train', 'Walk'))
ax.text(1, 1, "6")
plt.show()

<Figure size 288x288 with 0 Axes>

# 4.4(a): Binomial distribution.
from scipy.stats import binom
import matplotlib.pyplot as plt
n , p = 25 , 0.3
r_values = list(range(n+1))
dist = [binom.pmf(r , n , p) for r in r_values]
plt.bar(r_values , dist)
plt.show()
print("Mean = ", binom.mean(n = 25  , p = 0.3))
print("Variance = ", binom.var(n = 25 , p = 0.3))
print("Standard Deviation = ", binom.std(n = 25 , p = 0.3))

Mean =  7.5
Variance =  5.25
Standard Deviation =  2.29128784747792

# 4.4(b): Hypothesis testing.
n , p = 5 , 0.6
# H0: p = 0.6
# H1: p < 0.6
r_values = list(range(n+1))
dist = [binom.pmf(r , n , p) for r in r_values]
plt.bar(r_values , dist)
plt.show()
print("Mean = ", binom.mean(n = 5 , p = 0.6))
print("Variance = ", binom.var(n = 5 , p = 0.6))
print("Standard Deviation = ", binom.std(n = 5 , p = 0.6))
PXLTE1 = binom.cdf(1, n, p)  # Cumulative distribution function.
print("P(X<=1)=" , PXLTE1)      # P(X<=1) = 0.08704 > 0.05.
print("Do not reject H0. Insufficient evidence to suggest Thalia is correct.")

Mean =  3.0
Variance =  1.2000000000000002
Standard Deviation =  1.0954451150103324
P(X<=1)= 0.08704
Do not reject H0. Insufficient evidence to suggest Thalia is correct.

# 4.4(c): Velocity time graph.
xs = [0, 4, 6, 10 , 16]
ys = [2, 6, 2, 14, 18]
plt.axis([0, 18, 0, 20])
plt.plot(xs, ys)
plt.xlabel("Time(s)")
plt.ylabel("Velocity (m/s)")
plt.show()

# 4.4(d): Simultaneous equations.
a , b = symbols("a b")
solve([2 + 4 - 3 + a , 3 - 1 + 2 + b] , [a , b])

{a: -3, b: -4}

# 4.4(e): Velocity-time graph.
pts = np.array([0,2,3,4,4,4,6,1,3,5,7,7,1,0])
fig = plt.figure()
plt.plot(pts)
plt.xlabel("Time(s)")
plt.ylabel("Velocity (m/s)")
plt.show()

# 4.4(f): Proof.
# For the actual proof, use the principle of mathematical induction.
for n in range(0, 101):
  if ((2 * n + 3)**2 - (2 * n - 3)**2) % 6 is not 0:
    print("Failure")
  break
print("(2 * n + 3)**2 - (2 * n - 3)**2, is a multiple of 6, for 0<=n<=100.")

(2 * n + 3)**2 - (2 * n - 3)**2, is a multiple of 6, for 0<=n<=100.

Chapter 5: Python for A-Level Mathematics

# 5.1(a): Series expansion.
from sympy import *
x = symbols("x")
sec(x).series(x , 0 , 6)

$\displaystyle 1 + \frac{x^{2}}{2} + \frac{5 x^{4}}{24} + O\left(x^{6}\right)$

# 5.1(b): Summation of geometric progression.
n = symbols("n")
summation(2 * 3**(n - 1) , (n , 1, 50))

$\displaystyle 717897987691852588770248$

# 5.1(c): Functions of functions.
def f(x):
  return(2 * x**2 -1)
def g(x):
  return(x**3 + 1)
print("g(f(-1))=" , g(f(-1)))

g(f(-1))= 2

# 5.1(d): Differentiation.
x = symbols("x")
df = diff(2 * x / (x**2 + 3) , x)
df.subs(x , 2)

$\displaystyle - \frac{2}{49}$

# 5.1(e): Finding roots of trigonometric functions.
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0 , 2 * np.pi)
y1 = 3 * np.sin(t) + 4 * np.cos(t) - 5
y2 = 0 * t
plt.plot(t , y1 , t , y2 , 'r--' )
sol1 = nsolve(3 * sin(x) + 4 * cos(x) - 5 , (0 , 1))
print(sol1)

0.643501499704454

# 5.2(a): Partial fractions.
apart(3 * x / ((x+1)*(x-1)*(2*x-5)))

$\displaystyle \frac{10}{7 \left(2 x - 5\right)} - \frac{3}{14 \left(x + 1\right)} - \frac{1}{2 \left(x - 1\right)}$

# 5.2(b): Equality of functions.
t = np.linspace(0 , 2 * np.pi , 1000)
y1 = np.tan(2 * t)
y2 = 2 * np.tan(t) / (1 - np.tan(t)**2)
plt.ylim(-3,3)
plt.plot(t, y1 , "b-", linewidth = 4)
plt.plot(t , y2 , "r-", linewidth = 2)

[<matplotlib.lines.Line2D at 0x7fd3be2580a0>]

# 5.2(c): Differentiation.
diff(x * exp(sin(x)) / cos(x) , x)

$\displaystyle \frac{x e^{\sin{\left(x \right)}} \sin{\left(x \right)}}{\cos^{2}{\left(x \right)}} + x e^{\sin{\left(x \right)}} + \frac{e^{\sin{\left(x \right)}}}{\cos{\left(x \right)}}$

# 5.2(d): Simple integration.
print(integrate(x**3 * sin(x))," + c")

-x**3*cos(x) + 3*x**2*sin(x) + 6*x*cos(x) - 6*sin(x)  + c

# 5.2(e): Parametric plot.
t = np.linspace(0, 2 * np.pi, 100)
x = 4 * np.cos(t + np.pi / 6)
y = 2 * np.sin(t)
plt.plot(x, y)
plt.show()

# 5.3(a): Three-dimensional vectors.
u , v = np.array([-1 , 2 , 3]) , np.array([9 , 0 , -2]) 
np.linalg.norm(2 * u - 3 * v)

31.63858403911275

# 5.3(b): Differential equation, amount of drug in blood.
# Use dsolve and plot from the sympy library.
from sympy import *
t = symbols("t")
N = Function("N")
ODE = N(t).diff(t) + N(t) * (t + 1) / 8
Nsoln = dsolve(ODE, ics = {N(0):30})
print("N(t) = " , Nsoln.rhs)
plot(Nsoln.rhs , (t , 0 , 8), ylabel = "N(t)")
print("N(8) = " , Nsoln.subs(t , 8.0).rhs , "mg")

N(t) =  30*exp(t*(-t - 2)/16)

N(8) =  0.202138409972564 mg

# 5.3(c): Newton-Raphson method. Work to 8 significant figures throughout.
# You can determine dfn(x) using diff(x**2 * log(x) - 5 , x).
def fn(x):
  return x**2 * log(x) - 5
def dfn(x):
  return 2 * x * log(x) + x
def NewtonRaphson(x): 
  i=0
  h = round(fn(x) / dfn(x) , 8)
  while abs(h) >= 0.0001:
    h = round(fn(x) / dfn(x) , 8)
    x = x - h
    i += 1
    print("x(", i ,")=", x)
# Start at x( 0 ) = 2.
NewtonRaphson(2)

x( 1 )= 2.4667092
x( 2 )= 2.3953696
x( 3 )= 2.3935186
x( 4 )= 2.3935174

# 5.3(d): Probability density function (PDF) and probability.
k , x = symbols("k x")
sol = solve(64 * k - k / 4 - 1, k)
print("k=" , sol)
t = np.linspace(1 , 4 , 1000)
plt.plot(t , 4 * t**3 / 255 , label = "PDF")
plt.legend()
plt.show()
print("P(1<X<2)=",integrate(4 * x**3 / 255 , (x , 1 , 2 )))

k= [4/255]

P(1<X<2)= 1/17

# 5.3(e): Probability density function (PDF) and cumulative distribution function (CDF).
t = np.linspace(0 , 2 , 1000)
plt.plot(t , 3 * t**2 / 8 , label = "PDF")
plt.plot(t , t**3 / 8 , label = "CDF")
plt.legend()
plt.show()

# 5.3(f): Hypothesis testing.
# H0: mu = 80 and H1: mu > 80.
# Evidence: xbar=83 and n=100.
# Use the z-test: z=(xbar-mu) / (sigma/sqrt(n)).
# The critical value (-1.645) for a directional (upper-tail critical) 
# test at a 0.05 level of significance.
z = (83 - 80) / (15 / sqrt(100))
print("z=",z)
import scipy.stats as stats
import math
mu , sigma = 80 , 15
x = np.linspace(mu - 3*sigma, mu + 3*sigma, 100)
plt.plot(x, stats.norm.pdf(x, mu, sigma))
plt.show()
print("z > 1.645, reject H0, and the managers claim is supported.")

z= 2

z > 1.645, reject H0, and the managers claim is supported.

# 5.4(a): SUVAT
# v^2=u^2+2as.
print("Distance=" , (14**2-0**2) / (2 * 9.8),"m")
# v = u + at.
print("Time=" , (14 - 0) / 9.8 , "s")

Distance= 10.0 m
Time= 1.4285714285714284 s

# 5.4(b): Magnitude of resultant force.
FR = np.array([2-4+6 , -3+5+3])
print("Magnitude of FR=" , np.linalg.norm(FR) , "N")

Magnitude of FR= 6.4031242374328485 N

# 5.4(c): Forces in equilibrium.
R , T = symbols("R T")
Eq1 = R + T - 30*9.8 - 40*9.8
Eq2 = 2*40*9.8 + 3*30*9.8 - 5*R
print("Forces=",solve([Eq1 , Eq2], [R , T]), "N")

Forces= {R: 333.200000000000, T: 352.800000000000} N

# 5.4(d): Projectile clearing a wall.
# sx = vx0 * t and sy = vy0 * t - g * t^2 / 2.
t = 100 / 20
vy0 = (10 +9.8 * t**2 / 2) / t
print("Velocity in y-direction is", vy0 , "m/s")

Velocity in y-direction is 26.5 m/s

# 5.4(e): Mass on a plane.
# mu*m*g*cos(30)-m*g*sin(30)=0.
m , g = 2 , 9.8
mu = tan(np.radians(30))
print("mu=", mu)
# If the plane is inclined at 40 degrees.
# Parallel to the plane: m*g*sin(40) - mu * m * g * cos(40) = m * a.
a = (m * g *  np.sin(np.radians(40)) - mu * m * g * np.cos(np.radians(40))) / m
print("Acceleration down the slope =", a , "m/s/s")

mu= 0.577350269189626
Acceleration down the slope = 1.96501411355769 m/s/s

End of Solutions for Section 1