Python Doing Math

Python Doing Math

written by sean base on following books

Github | https://github.com/newsteinking/High_DoingMath

Chapter 0: About

Python Advanced Curriculm

by sean

Base on Python Doing Math

Thanks to

잠시나마 오픈소스에 대한 희망을 안겨주었던 멤버들에게 감사의 마음을 전합니다.
다들 다른곳에서 각자 열심히 일하는데 모두들 건승하길 바랍니다.

.

• sean
• Mr Ju SS
• OSS Members

SEAN’s Paradise

I think that My Life as Software Engineer was torrible , but it’s role for social is important so, I keep going for better life & software development

chapter 1: Working with numbers

1.1 Basic Mathematical Operations

… code-block:: python

>>> 1 + 2
3
>>> 1 + 3.5
4.5
>>> -1 + 2.5
1.5
>>> 100 – 45
55
>>> -1.1 + 5
3.9

>>> 3 * 2
6
>>> 3.5 * 1.5
5.25

>>> 3 / 2
1.5
>>> 4 / 2
2.0

>>> 3 // 2
1

>>> -3 // 2
-2

>>> 9 % 2
1

>>> 2 ** 2
4
>>> 2 ** 10
1024
>>> 1 ** 10
1

>>> 8 ** (1/3)
2.0

>>> 5 + 5 * 5
30
>>> (5 + 5) * 5
50

1.2 Labels: Attaching Names to Numbers

… code-block:: python

>>> a = 3
>>> a + 1
4
>>> a = 5
>>> a + 1
6

1.3 Different Kinds of Numbers

… code-block:: python

>>> type(3)
<class 'int'>
>>> type(3.5)
<class 'float'>
>>> type(3.0)
<class 'float'>

>>> int(3.8)
3
>>> int(3.0)
3

>>> float(3)
3.0

1.4 Working with Fractions

… code-block:: python

u >>> from fractions import Fraction v >>> f = Fraction(3, 4) w >>> f Fraction(3, 4)

>>> Fraction(3, 4) + 1 + 1.5
3.25

>>> Fraction(3, 4) + 1 + Fraction(1/4)
Fraction(2, 1)

1.5 Complex Numbers

… code-block:: python

>>> a = 2 + 3j
>>> type(a)
<class 'complex'>

>>> a = complex(2, 3)
>>> a
(2 + 3j)

>>> b = 3 + 3j
>>> a + b
(5 + 6j)
>>> a - b
(-1 + 0j)

>>> a * b
(-3 + 15j)
>>> a / b
(0.8333333333333334 + 0.16666666666666666j)

>>> z = 2 + 3j
>>> z.real
2.0
>>> z.imag
3.0

>>> z.conjugate()
(2 - 3j)

>>> (z.real ** 2 + z.imag ** 2) ** 0.5
3.605551275463989

>>> abs(z)
3.605551275463989

1.6 Getting User Input

… code-block:: python

u >>> a = input() v 1

>>> a
w '1'

>>> s1 = 'a string'
>>> s2 = "a string"

>>> a = '1'
>>> int(a) + 1
2
>>> float(a) + 1
2.0

>>> int('2.0')
Traceback (most recent call last):
File "<pyshell#26>", line 1, in <module>
int('2.0')
ValueError: invalid literal for int() with base 10: '2.0'

>>> a = float(input())
3/4
Traceback (most recent call last):
File "<pyshell#25>", line 1, in <module>
a=float(input())
ValueError: could not convert string to float: '3/4'

1.8 Handling Exceptions and Invalid Input

… code-block:: python

>>> try:
a = float(input('Enter a number: '))
except ValueError:
print('You entered an invalid number')

Enter a number: 3/4 u You entered an invalid number

>>> a = input('Input an integer: ')

>>> a = int(input())
1
>>> a + 1
2

>>> a = int(input())
1.0
Traceback (most recent call last):
File "<pyshell#42>", line 1, in <module>
a=int(input())
ValueError: invalid literal for int() with base 10: '1.0'

>>> 1.1.is_integer()
False

>>> 1.0.is_integer()
True

1.9 Fractions and Complex Numbers as Input

… code-block:: python

>>> a = Fraction(input('Enter a fraction: '))
Enter a fraction: 3/4
>>> a
Fraction(3, 4)

>>> a = Fraction(input('Enter a fraction: '))
Enter a fraction: 3/0
Traceback (most recent call last):
File "<pyshell#2>", line 1, in <module>
a = Fraction(input('Enter a fraction: '))
File "/usr/lib64/python3.3/fractions.py", line 167, in __new__
raise ZeroDivisionError('Fraction(%s, 0)' % numerator)
ZeroDivisionError: Fraction(3, 0)

>>> try:
a = Fraction(input('Enter a fraction: '))
except ZeroDivisionError:
print('Invalid fraction')
Enter a fraction: 3/0
Invalid fraction

>>> z = complex(input('Enter a complex number: '))
Enter a complex number: 2+3j
>>> z
(2+3j)

>>> z = complex(input('Enter a complex number: '))
Enter a complex number: 2 + 3j
Traceback (most recent call last):
File "<pyshell#43>", line 1, in <module>
z = complex(input('Enter a complex number: '))
ValueError: complex() arg is a malformed string

1.10 Writing Programs That Do the Math for You

Calculating the Factors of an Integer

… code-block:: python

>>> def is_factor(a, b):
if b % a == 0:
return True
else:
return False

>>> is_factor(4, 1024)
True

>>> for i in range(1, 4):
print(i)

>>> for i in range(5):
print(i)

>>> for i in range(1,10,2):
print(i)

‘’’ Find the factors of an integer ‘’’ def factors(b): u for i in range(1, b+1): if b % i == 0: print(i) if __name__ == ‘__main__’: b = input(‘Your Number Please: ‘) b = float(b) v if b > 0 and b.is_integer(): factors(int(b)) else: print(‘Please enter a positive integer’)

Your Number Please: 25 1 5 25

Your Number Please: 15.5 Please enter a positive integer

Generating Multiplication Tables

… code-block:: python

‘’’ enhanced_multi_table.py

Multiplication table printer: Enter the number and the number of multiples to be printed ‘’‘

def multi_table(a, n):

for i in range(1, n+1):

print(‘{0} x {1} = {2}’.format(a, i, a*i))

if __name__ == ‘__main__’:

try:

a = float(input(‘Enter a number: ‘)) n = float(input(‘Enter the number of multiples: ‘)) if not n.is_integer() or n < 0:

print(‘The number of multiples should be a positive integer’)

else:

multi_table(a, int(n))

except ValueError:

print(‘You entered an invalid input’)

Converting Units of Measurement

… code-block:: python

‘’’ enhanced_unit_converter_exit_power.py

Unit converter:

Miles and Kilometers Kilograms and Pounds Celsius and Fahrenheit

‘’‘

def print_menu():

print(‘1. Kilometers to Miles’) print(‘2. Miles to Kilometers’) print(‘3. Kilograms to Pounds’) print(‘4. Pounds to Kilograms’) print(‘5. Celsius to Fahrenheit’) print(‘6. Fahrenheit to Celsius’)

def km_miles():

km = float(input(‘Enter distance in kilometers: ‘)) miles = km / 1.609

print(‘Distance in miles: {0}’.format(miles))

def miles_km():

miles = float(input(‘Enter distance in miles: ‘)) km = miles * 1.609

print(‘Distance in kilometers: {0}’.format(km))

def kg_pounds():

kg = float(input(‘Enter weight in kilograms: ‘)) pounds = kg * 2.205

print(‘Weight in pounds: {0}’.format(pounds))

def pounds_kg():

pounds = float(input(‘Enter weight in pounds: ‘)) kg = pounds / 2.205

print(‘Weight in kilograms: {0}’.format(kg))

def cel_fahren():

celsius = float(input(‘Enter temperature in celsius: ‘)) fahrenheit = celsius* (9 / 5) + 32 print(‘Temperature in fahrenheit: {0}’.format(fahrenheit))

def fahren_cel():

fahrenheit = float(input(‘Enter temperature in fahrenheit: ‘)) celsius = (fahrenheit - 32)*(5/9) print(‘Temperature in celsius: {0}’.format(celsius))

if __name__ == ‘__main__’:

print_menu()

while True:

choice = input(‘Which conversion would you like to do? ‘)

if choice == ‘1’:

km_miles()

if choice == ‘2’:

miles_km()

if choice == ‘3’:

kg_pounds()

if choice == ‘4’:

pounds_kg()

if choice == ‘5’:

cel_fahren()

if choice == ‘6’:

fahren_cel()

answer = input(‘Do you want to exit? (y) for yes ‘) if answer == ‘y’:

break

Finding the Roots of a Quadratic Equation

… code-block:: python

‘’’ Quadratic equation root calculator ‘’’ def roots(a, b, c):

D = (b*b - 4*a*c)**0.5 x_1 = (-b + D)/(2*a) x_2 = (-b - D)/(2*a) print(‘x1: {0}’.format(x_1)) print(‘x2: {0}’.format(x_2))

if __name__ == ‘__main__’:

a = input(‘Enter a: ‘) b = input(‘Enter b: ‘) c = input(‘Enter c: ‘) roots(float(a), float(b), float(c))

chapter 2: visualizing Data with Graphs

이 장에서는 다양한 데이터를 mathplotlib를 통해 표현하는 방법을 배우도록 하겠다.

2.1 Understanding the Cartesian Coordinate Plane

다음처럼 숫자가 새겨진 라인을 생각해 보자.

아래는 Cartesian coordinate plane 을 나타낸다.

Working with Lists and Tuples

그래프를 그릴때 list,tuple을 많이 사용하게 될것이다.

Iterating over a List or Tuple

>>> l = [1, 2, 3]
>>> for item in l:
print(item)

>>> l = [1, 2, 3]
>>> for index, item in enumerate(l):
print(index, item)

Creating Graphs with Matplotlib

>>> x_numbers = [1, 2, 3]
>>> y_numbers = [2, 4, 6]

>>> from pylab import plot, show
>>> plot(x_numbers, y_numbers)
[<matplotlib.lines.Line2D object at 0x7f83ac60df10>]

plot(x_numbers, y_numbers, marker='o')

plot(x_numbers, y_numbers, 'o')

Graphing the Average Annual Temperature in New York City

from pylab import plot, show

nyc_temp = [53.9, 56.3, 56.4, 53.4, 54.5, 55.8, 56.8, 55.0, 55.3, 54.0, 56.7, 56.4, 57.3]


years = range(2000, 2013)



#plot(nyc_temp, marker='o')

years = range(2000, 2013)
plot(years, nyc_temp, marker='o')

show()

Comparing the Monthly Temperature Trends of New York City

from pylab import plot,show
from pylab import legend
legend([2000, 2006, 2012])


nyc_temp_2000 = [31.3, 37.3, 47.2, 51.0, 63.5, 71.3, 72.3, 72.7, 66.0, 57.0, 45.3, 31.1]
nyc_temp_2006 = [40.9, 35.7, 43.1, 55.7, 63.1, 71.0, 77.9, 75.8, 66.6, 56.2, 51.9, 43.6]
nyc_temp_2012 = [37.3, 40.9, 50.9, 54.8, 65.1, 71.0, 78.8, 76.7, 68.8, 58.0, 43.9, 41.5]

months = range(1, 13)

#plot(months, nyc_temp_2000, months, nyc_temp_2006, months, nyc_temp_2012)

#===============================================================================
# plot(months, nyc_temp_2000)
# plot(months, nyc_temp_2006)
# plot(months, nyc_temp_2012)
#===============================================================================

plot(months, nyc_temp_2000, months, nyc_temp_2006, months, nyc_temp_2012)

show()

Customizing Graphs

Adding a Title and Labels

from pylab import plot, show, title, xlabel, ylabel, legend

months = range(1, 13)

nyc_temp_2000 = [31.3, 37.3, 47.2, 51.0, 63.5, 71.3, 72.3, 72.7, 66.0, 57.0, 45.3, 31.1]
nyc_temp_2006 = [40.9, 35.7, 43.1, 55.7, 63.1, 71.0, 77.9, 75.8, 66.6, 56.2, 51.9, 43.6]
nyc_temp_2012 = [37.3, 40.9, 50.9, 54.8, 65.1, 71.0, 78.8, 76.7, 68.8, 58.0, 43.9, 41.5]


plot(months, nyc_temp_2000, months, nyc_temp_2006, months, nyc_temp_2012)

title('Average monthly temperature in NYC')
xlabel('Month')
ylabel('Temperature')
legend([2000, 2006, 2012])

show()

Customizing the Axes

from pylab import plot, show, axis

nyc_temp = [53.9, 56.3, 56.4, 53.4, 54.5, 55.8, 56.8, 55.0, 55.3, 54.0, 56.7, 56.4, 57.3]
plot(nyc_temp, marker='o')
print(axis())
axis(ymin=0)

show()

Plotting Using pyplot

   '''
   Simple plot using pyplot
   '''
   import matplotlib.pyplot
   def create_graph():
       x_numbers = [1, 2, 3]
       y_numbers = [2, 4, 6]
       matplotlib.pyplot.plot(x_numbers, y_numbers)
       matplotlib.pyplot.show()

   if __name__ == '__main__':
       create_graph()

Saving the Plots

---

from pylab import plot, savefig
x = [1, 2, 3]
y = [2, 4, 6]
plot(x, y)
savefig('mygraph.png')

#savefig('C:\mygraph.png')

Plotting with Formulas

Newton’s Law of Universal Gravitation

Newton’s law of universal gravitation

'''
The relationship between gravitational force and
distance between two bodies
'''
import matplotlib.pyplot as plt
# Draw the graph
def draw_graph(x, y):
    plt.plot(x, y, marker='o')
    plt.xlabel('Distance in meters')
    plt.ylabel('Gravitational force in newtons')
    plt.title('Gravitational force and distance')
    plt.show()

def generate_F_r():
# Generate values for r
    r = range(100, 1001, 50)
# Empty list to store the calculated values of F
    F = []
# Constant, G
    G = 6.674*(10**-11)
# Two masses
    m1 = 0.5
    m2 = 1.5
# Calculate force and add it to the list, F
    for dist in r:
        force = G*(m1*m2)/(dist**2)
        F.append(force)
# Call the draw_graph function
    draw_graph(r, F)
if __name__=='__main__':
    generate_F_r()

Projectile Motion

Generating Equally Spaced Floating Point Numbers

'''
Generate equally spaced floating point
numbers between two given values
'''
def frange(start, final, increment):
    numbers = []
    while start < final:
        numbers.append(start)
        start = start + increment
        return numbers

Drawing the Trajectory

'''
Draw the trajectory of a body in projectile motion
'''
from matplotlib import pyplot as plt
import math
def draw_graph(x, y):
    plt.plot(x, y)
    plt.xlabel('x-coordinate')
    plt.ylabel('y-coordinate')
    plt.title('Projectile motion of a ball')
def frange(start, final, interval):
    numbers = []
    while start < final:
        numbers.append(start)
        start = start + interval
        return numbers
def draw_trajectory(u, theta):
    theta = math.radians(theta)
    g = 9.8
# Time of flight
    t_flight = 2*u*math.sin(theta)/g
# Find time intervals
    intervals = frange(0, t_flight, 0.001)

    # List of x and y coordinates
    x = []
    y = []
    for t in intervals:
        x.append(u*math.cos(theta)*t)
        y.append(u*math.sin(theta)*t - 0.5*g*t*t)
        draw_graph(x, y)
if __name__ == '__main__':
    try:
        u = float(input('Enter the initial velocity (m/s): '))
        theta = float(input('Enter the angle of projection (degrees): '))
    except ValueError:
        print('You entered an invalid input')
    else:
        draw_trajectory(u, theta)
        plt.show()

Comparing the Trajectory at Different Initial Velocities

'''
Draw the trajectory of a body in projectile motion
'''
from matplotlib import pyplot as plt
import math
def draw_graph(x, y):
    plt.plot(x, y)
    plt.xlabel('x-coordinate')
    plt.ylabel('y-coordinate')
    plt.title('Projectile motion of a ball')
def frange(start, final, interval):
    numbers = []
    while start < final:
        numbers.append(start)
        start = start + interval
        return numbers
def draw_trajectory(u, theta):
    theta = math.radians(theta)
    g = 9.8
# Time of flight
    t_flight = 2*u*math.sin(theta)/g
# Find time intervals
    intervals = frange(0, t_flight, 0.001)

    # List of x and y coordinates
    x = []
    y = []
    for t in intervals:
        x.append(u*math.cos(theta)*t)
        y.append(u*math.sin(theta)*t - 0.5*g*t*t)
        draw_graph(x, y)
if __name__ == '__main__':
    # List of three different initial velocities
    u_list = [20, 40, 60]
    theta = 45
    for u in u_list:
        draw_trajectory(u, theta)
# Add a legend and show the graph
    plt.legend(['20', '40', '60'])
    plt.show()

chapter 3: Describing Data With Statstics

이 장에서는 통계치를 가지고 한번 파이썬 코드를 짜보자.

3.1 Finding the Mean

>>> shortlist = [1, 2, 3]
>>> sum(shortlist)
6

>>> len(shortlist)
3

'''
Calculating the mean
'''
def calculate_mean(numbers):
    s = sum(numbers)
    N = len(numbers)
# Calculate the mean
    mean = s/N
    return mean


if __name__ == '__main__':
    donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
    mean = calculate_mean(donations)
    N = len(donations)
    print('Mean donation over the last {0} days is {1}'.format(N, mean))

Finding the Median

samplelist = [4, 1, 3]
samplelist.sort()
samplelist
[1, 3, 4]

'''
Calculating the median
'''
def calculate_median(numbers):
    N = len(numbers)
    numbers.sort()
    # Find the median
    if N % 2 == 0:
        # if N is even
        m1 = N/2
        m2 = (N/2) + 1
        # Convert to integer, match position
        m1 = int(m1) - 1
        m2 = int(m2) - 1
        median = (numbers[m1] + numbers[m2])/2
    else:
        m = (N+1)/2
        # Convert to integer, match position
        m = int(m) - 1
        median = numbers[m]
    return median

if __name__ == '__main__':
    donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
    median = calculate_median(donations)
    N = len(donations)
    print('Median donation over the last {0} days is {1}'.format(N, median))

Finding the Mode and Creating a Frequency Table

Finding the Most Common Elements

>>> simplelist = [4, 2, 1, 3, 4]
>>> from collections import Counter
>>> c = Counter(simplelist)
>>> c.most_common()
[(4, 2), (1, 1), (2, 1), (3, 1)]

>>> c.most_common(1)
[(4, 2)]

>>> c.most_common(2)
[(4, 2), (1, 1)]

>>> mode = c.most_common(1)
>>> mode
[(4, 2)]
v >>> mode[0]
(4, 2)
w >>> mode[0][0]
4

Finding the Mode

'''
Calculating the mode
'''
from collections import Counter
def calculate_mode(numbers):
    c = Counter(numbers)
    mode = c.most_common(1)
    return mode[0][0]
if __name__=='__main__':
    scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
    mode = calculate_mode(scores)
    print('The mode of the list of numbers is: {0}'.format(mode))


'''
Calculating the mode when the list of numbers may
have multiple modes
'''
from collections import Counter
def calculate_mode(numbers):
    c = Counter(numbers)
    numbers_freq = c.most_common()
    max_count = numbers_freq[0][1]
    modes = []
    for num in numbers_freq:
        if num[1] == max_count:
            modes.append(num[0])
    return modes
if __name__ == '__main__':
    scores = [5, 5, 5, 4, 4, 4, 9, 1, 3]
    modes = calculate_mode(scores)
    print('The mode(s) of the list of numbers are:')
    for mode in modes:
        print(mode)

Creating a Frequency Table

'''
Frequency table for a list of numbers
'''
from collections import Counter

def frequency_table(numbers):
    table = Counter(numbers)
    print('Number\tFrequency')
    for number in table.most_common():
        print('{0}\t{1}'.format(number[0], number[1]))

if __name__=='__main__':
    scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
    frequency_table(scores)



'''
Frequency table for a list of numbers
Enhanced to display the table sorted by the numbers
'''
from collections import Counter
def frequency_table(numbers):
    table = Counter(numbers)
    numbers_freq = table.most_common()
    numbers_freq.sort()
    print('Number\tFrequency')
    for number in numbers_freq:
        print('{0}\t{1}'.format(number[0], number[1]))

if __name__ == '__main__':
    scores = [7, 8, 9, 2, 10, 9, 9, 9, 9, 4, 5, 6, 1, 5, 6, 7, 8, 6, 1, 10]
    frequency_table(scores)

Measuring the Dispersion

Finding the Range of a Set of Numbers

'''
Find the range
'''
def find_range(numbers):
    lowest = min(numbers)
    highest = max(numbers)
    # Find the range
    r = highest-lowest
    return lowest, highest, r
if __name__ == '__main__':
    donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
    lowest, highest, r = find_range(donations)
    print('Lowest: {0} Highest: {1} Range: {2}'.format(lowest, highest, r))

Finding the Variance and Standard Deviation

'''
Find the variance and standard deviation of a list of numbers
'''
def calculate_mean(numbers):
    s = sum(numbers)
    N = len(numbers)
    # Calculate the mean
    mean = s/N
    return mean
def find_differences(numbers):
    # Find the mean
    mean = calculate_mean(numbers)
    # Find the differences from the mean
    diff = []

    for num in numbers:
        diff.append(num-mean)
    return diff

    def calculate_variance(numbers):
        # Find the list of differences
        diff = find_differences(numbers)
        # Find the squared differences
        squared_diff = []
        for d in diff:
            squared_diff.append(d**2)
        # Find the variance
        sum_squared_diff = sum(squared_diff)
        variance = sum_squared_diff/len(numbers)
        return variance

if __name__ == '__main__':
    donations = [100, 60, 70, 900, 100, 200, 500, 500, 503, 600, 1000, 1200]
    variance = calculate_variance(donations)
    print('The variance of the list of numbers is {0}'.format(variance))
    std = variance**0.5
    print('The standard deviation of the list of numbers is {0}'.format(std))

Calculating the Correlation Between Two Data Sets

def find_corr_x_y(x,y):
    n = len(x)
    # Find the sum of the products
    prod = []
    for xi,yi in zip(x,y):
        prod.append(xi*yi)
    sum_prod_x_y = sum(prod)
    sum_x = sum(x)
    sum_y = sum(y)
    squared_sum_x = sum_x**2
    squared_sum_y = sum_y**2
    x_square = []
    for xi in x:
        x_square.append(xi**2)
    # Find the sum
    x_square_sum = sum(x_square)
    y_square=[]
    for yi in y:
        y_square.append(yi**2)
    # Find the sum
    y_square_sum = sum(y_square)

    # Use formula to calculate correlation
    numerator = n*sum_prod_x_y - sum_x*sum_y
    denominator_term1 = n*x_square_sum - squared_sum_x
    denominator_term2 = n*y_square_sum - squared_sum_y
    denominator = (denominator_term1*denominator_term2)**0.5
    correlation = numerator/denominator
    return correlation

High School Grades and Performance on College Admission Tests

Reading Data from Files

# Find the sum of numbers stored in a file
def sum_data(filename):
    s = 0
    with open(filename) as f:
        for line in f:
            s = s + float(line)
        print('Sum of the numbers: {0}'.format(s))
if __name__ == '__main__':
    sum_data('mydata.txt')


'''
Calculating the mean of numbers stored in a file
'''
def read_data(filename):
    numbers = []
    with open(filename) as f:
        for line in f:
            numbers.append(float(line))
    return numbers
def calculate_mean(numbers):
    s = sum(numbers)
    N = len(numbers)
    mean = s/N
    return mean
if __name__ == '__main__':
    data = read_data('mydata.txt')
    mean = calculate_mean(data)
    print('Mean: {0}'.format(mean))

Reading Data from a CSV File

import csv
import matplotlib.pyplot as plt
def scatter_plot(x, y):
    plt.scatter(x, y)
    plt.xlabel('Number')
    plt.ylabel('Square')
    plt.show()
def read_csv(filename):
    numbers = []
    squared = []
    with open(filename) as f:
        reader = csv.reader(f)
        next(reader)
        for row in reader:
            numbers.append(int(row[0]))
            squared.append(int(row[1]))
        return numbers, squared
if __name__ == '__main__':
    numbers, squared = read_csv('numbers.csv')
    scatter_plot(numbers, squared)

chapter 4: Algebra And Symbolic Math with Sympy

이 장에서는 SymPy 라이브러리를 집중 분석해 보도록 하자.

4.1 Defining Symbols and Symbolic Operations

>>> from sympy import Symbol
>>> x = Symbol('x')

>>> from sympy import Symbol
>>> x = Symbol('x')
>>> x + x + 1
2*x + 1


>>> x = Symbol('x')
>>> y = Symbol('y')
>>> z = Symbol('z')


>>> from sympy import symbols
>>> x,y,z = symbols('x,y,z')

>>> from sympy import Symbol
>>> x = Symbol('x')
>>> y = Symbol('y')
96 Chapter 4
>>> s = x*y + x*y
>>> s
2*x*y

>>> p = x*(x + x)
>>> p
2*x**2

>>> p = (x + 2)*(x + 3)
>>> p
(x + 2)*(x + 3)

4.2 Working with Expressions

Factorizing and Expanding Expressions

x2 − y2 = (x + y)(x − y).

from sympy import Symbol
x = Symbol('x')
y = Symbol('y')

from sympy import factor,expand
expr = x**2 - y**2
print(factor(expr))
#(x - y)*(x + y)

factors = factor(expr)
print(expand(factors))
#x**2 - y**2

x3 + 3x2y + 3xy2 + y3 = (x + y)3

from sympy import factor,expand
expr = x**3 + 3*x**2*y + 3*x*y**2 + y**3
factors = factor(expr)
print(factors)
#(x + y)**3
print(expand(factors))
#x**3 + 3*x**2*y + 3*x*y**2 + y**3

Pretty Printing

from sympy import Symbol,factor,expand
x=Symbol('x')
y=Symbol('y')
expr = x*x + 2*x*y + y*y
print(expr)


from sympy import Symbol,factor,expand
x=Symbol('x')
y=Symbol('y')
expr = x*x + 2*x*y + y*y
print(expr)

from sympy import pprint

pprint(expr)

expr = 1 + 2*x + 2*x**2

pprint(expr)

#-*-coding:utf-8
from sympy import Symbol,factor,expand
x=Symbol('x')
y=Symbol('y')
expr = x*x + 2*x*y + y*y
print(expr)

from sympy import pprint

pprint(expr)

expr = 1 + 2*x + 2*x**2

pprint(expr)


from sympy import init_printing
init_printing(order='rev-lex')
pprint(expr)
#1 + 2·x + 2·x 2

Printing a Series

'''
Print the series:
x + x**2 + x**3 + ... + x**n
____ _____ _____
2 3 n
'''
from sympy import Symbol, pprint, init_printing
def print_series(n):
    # Initialize printing system with reverse order
    init_printing(order='rev-lex')
    x = Symbol('x')
    series = x
    for i in range(2, n+1):
        series = series + (x**i)/i
    pprint(series)

if __name__ == '__main__':
        n = input('Enter the number of terms you want in the series: ')
        print_series(int(n))

Substituting in Values

x2 + 2xy + y2

from sympy import Symbol

x=Symbol('x')
y=Symbol('y')

print(x*x + x*y + x*y + y*y)

expr = x*x + x*y + x*y + y*y

res = expr.subs({x:1, y:2})

print(res)

print(expr.subs({x:1-y}))


expr_subs = expr.subs({x:1-y})

from sympy import simplify

print(simplify(expr_subs))

Calculating the Value of a Series

'''
Print the series:
x + x**2 + x**3 + ... + x**n
____ _____ _____
2 3 n
'''
from sympy import Symbol, pprint, init_printing
def print_series(n, x_value):
    # Initialize printing system with reverse order
    init_printing(order='rev-lex')
    x = Symbol('x')
    series = x
    for i in range(2, n+1):
        series = series + (x**i)/i
    pprint(series)
    # Evaluate the series at x_value
    series_value = series.subs({x:x_value})
    print('Value of the series at {0}: {1}'.format(x_value, series_value))
if __name__ == '__main__':
    n = input('Enter the number of terms you want in the series: ')
    x_value = input('Enter the value of x at which you want to evaluate the series: ')
    print_series(int(n), float(x_value))

Converting Strings to Mathematical Expressions

from sympy import sympify
expr = input('Enter a mathematical expression: ')
#Enter a mathematical expression: x**2 + 3*x + x**3 + 2*x
expr = sympify(expr)

print(2*expr)

#expr = input('Enter a mathematical expression: ')
#Enter a mathematical expression: x**2 + 3*x + x**3 + 2x
#expr = sympify(expr)
#print(expr)


from sympy import sympify
from sympy.core.sympify import SympifyError
expr = input('Enter a mathematical expression: ')
#Enter a mathematical expression: x**2 + 3*x + x**3 + 2x
try:
    expr = sympify(expr)
except SympifyError:
    print('Invalid input')

Expression Multiplier

'''
Product of two expressions
'''
from sympy import expand, sympify
from sympy.core.sympify import SympifyError
def product(expr1, expr2):
    prod = expand(expr1*expr2)
    print(prod)
if __name__=='__main__':
    expr1 = input('Enter the first expression: ')
    #Enter the first expression: x**2 + x*2 + x
    expr2 = input('Enter the second expression: ')
    #Enter the second expression: x**3 + x*3 + x
    try:
        expr1 = sympify(expr1)
        expr2 = sympify(expr2)
    except SympifyError:
        print('Invalid input')
    else:
        product(expr1, expr2)

#x**5 + 3*x**4 + 4*x**3 + 12*x**2


#===============================================================================
# Enter the first expression: x*y+x
# Enter the second expression: x*x+y
# x**3*y + x**3 + x*y**2 + x*y
#===============================================================================

Solving Equations

from sympy import Symbol, solve
x = Symbol('x')
expr = x - 5 - 7
print(solve(expr))

Solving Quadratic Equations

from sympy import solve,Symbol
x = Symbol('x')
expr = x**2 + 5*x + 4
print(solve(expr, dict=True))
#x [{x: -4}, {x: -1}]

x2 + x + 1 = 0

   from sympy import Symbol,solve
   x=Symbol('x')
   expr = x**2 + x + 1
   print(solve(expr, dict=True))
   #[{x: -1/2 - sqrt(3)*I/2}, {x: -1/2 + sqrt(3)*I/2}]


Solving for One Variable in Terms of Others

---

ax2 + bx + c = 0

from sympy import symbol

x = Symbol('x')
a = Symbol('a')
b = Symbol('b')
c = Symbol('c')

expr = a*x*x + b*x + c
print(solve(expr, x, dict=True))
#[{x: (-b + sqrt(-4*a*c + b**2))/(2*a)}, {x: -(b + sqrt(-4*a*c + b**2))/(2*a)}]

   from sympy import Symbol, solve, pprint
   s = Symbol('s')
   u = Symbol('u')
   t = Symbol('t')
   a = Symbol('a')
   expr = u*t + (1/2)*a*t*t - s
   t_expr = solve(expr,t, dict=True)
   pprint(t_expr)


.. image:: ./img/chapter4-3.png

Solving a System of Linear Equations

2x + 3y = 6 3x + 2y = 12

#-*-coding=utf-8
from sympy import Symbol,solve

x = Symbol('x')
y = Symbol('y')
expr1 = 2*x + 3*y - 6
expr2 = 3*x + 2*y - 12

print(solve((expr1, expr2), dict=True))


soln = solve((expr1, expr2), dict=True)
soln = soln[0]
print(expr1.subs({x:soln[x], y:soln[y]}))
#0
print(expr2.subs({x:soln[x], y:soln[y]}))
#0

Plotting Using SymPy

from sympy.plotting import plot
from sympy import Symbol
x = Symbol('x')
#plot(2*x+3)

#plot((2*x + 3), (x, -5, 5))

#plot(2*x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3')

p = plot(2*x + 3, (x, -5, 5), title='A Line', xlabel='x', ylabel='2x+3', show=False)
p.save('line.png')

Plotting Expressions Input by the User

y = 2x + 3,

from sympy import Symbol,solve,sympify
expr = input('Enter an expression: ')
#Enter an expression: 2*x + 3*y - 6
expr = sympify(expr)
y = Symbol('y')
print(solve(expr, y))

'''
Plot the graph of an input expression
'''
from sympy import Symbol, sympify, solve
from sympy.plotting import plot
def plot_expression(expr):
    y = Symbol('y')
    solutions = solve(expr, y)
    expr_y = solutions[0]
    plot(expr_y)

if __name__=='__main__':
    expr = input('Enter your expression in terms of x and y: ')
    try:
        expr = sympify(expr)
    except SympifyError:
        print('Invalid input')
    else:
        plot_expression(expr)

Plotting Multiple Functions

from sympy.plotting import plot
from sympy import Symbol
x = Symbol('x')
plot(2*x+3, 3*x+1)

from sympy.plotting import plot
from sympy import Symbol
x = Symbol('x')
p = plot(2*x+3, 3*x+1, legend=True, show=False)
p[0].line_color = 'b'
p[1].line_color = 'r'
p.show()

chapter 5: Playing with Sets and Probability

이 장에서는 Probability 와 Random에 대해서 집중 공부하도록 하자.

5.1 What’s a Set?

Set Construction

Sympy를 이용하여 Python을 설정하는것을 알아보자.

from sympy import FiniteSet
s = FiniteSet(2, 4, 6)
print(s)


from sympy import FiniteSet
from fractions import Fraction
s = FiniteSet(1, 1.5, Fraction(1, 5))
print(s)

print(len(s))

print(4 in s)

s = FiniteSet()
print(s)
#EmptySet()


members = [1, 2, 3]
s = FiniteSet(*members)
print(s)
#{1, 2, 3}

Set Repetition and Order

from sympy import FiniteSet
members = [1, 2, 3, 2]
print(FiniteSet(*members))


from sympy import FiniteSet
s = FiniteSet(1, 2, 3)
for member in s:
    print(member)


from sympy import FiniteSet
s = FiniteSet(3, 4, 5)
t = FiniteSet(5, 4, 3)
print(s == t)

Subsets, Supersets, and Power Sets

from sympy import FiniteSet
s = FiniteSet(1)
t = FiniteSet(1,2)
print(s.is_subset(t))
#True
print(t.is_subset(s))
#False


print(s.is_superset(t))
#False
print(t.is_superset(s))
#True


s = FiniteSet(1, 2, 3)
ps = s.powerset()
print(ps)


from sympy import FiniteSet
s = FiniteSet(1, 2, 3)
t = FiniteSet(1, 2, 3)
print(s.is_proper_subset(t))
#False
print(t.is_proper_superset(s))
#False

t = FiniteSet(1, 2, 3, 4)
print(s.is_proper_subset(t))
#True
print(t.is_proper_superset(s))
#True

Set Operations

Union and Intersection

from sympy import FiniteSet
s = FiniteSet(1, 2, 3)
t = FiniteSet(2, 4, 6)
print(s.union(t))
#{1, 2, 3, 4, 6}


s = FiniteSet(1, 2)
t = FiniteSet(2, 3)
print(s.intersect(t))
#{2}


from sympy import FiniteSet
s = FiniteSet(1, 2, 3)
t = FiniteSet(2, 4, 6)
u = FiniteSet(3, 5, 7)
print(s.union(t).union(u))
#{1, 2, 3, 4, 5, 6, 7}


print(s.intersect(t).intersect(u))
#EmptySet()

Cartesian Product

from sympy import FiniteSet
s = FiniteSet(1, 2)
t = FiniteSet(3, 4)
p = s*t
print(p)

for elem in p:
    print(elem)


print(len(p) == len(s)*len(t))


from sympy import FiniteSet
s = FiniteSet(1, 2)
p = s**3
print(p)

for elem in p:
    print(elem)

Applying a Formula to Multiple Sets of Variables

from sympy import FiniteSet, pi
def time_period(length):
    g = 9.8
    T = 2*pi*(length/g)**0.5
    return T
if __name__ == '__main__':
    L = FiniteSet(15, 18, 21, 22.5, 25)
    for l in L:
        t = time_period(l/100)
        print('Length: {0} cm Time Period: {1:.3f} s'. format(float(l), float(t)))

Different Gravity, Different Results

from sympy import FiniteSet, pi
def time_period(length, g):
    T = 2*pi*(length/g)**0.5
    return T
if __name__ == '__main__':
    L = FiniteSet(15, 18, 21, 22.5, 25)
    g_values = FiniteSet(9.8, 9.78, 9.83)
    print('{0:^15}{1:^15}{2:^15}'.format('Length(cm)', 'Gravity(m/s^2)', 'Time Period(s)'))
    for elem in L*g_values:
        l = elem[0]
        g = elem[1]
        t = time_period(l/100, g)
        print('{0:^15}{1:^15}{2:^15.3f}'.format(float(l), float(g), float(t)))

Probability

def probability(space, event):
    return len(event)/len(space)
    def check_prime(number):
        if number != 1:
            for factor in range(2, number):
                if number % factor == 0:
                    return False
        else:
            return False
        return True
if __name__ == '__main__':
    space = FiniteSet(*range(1, 21))
    primes = []
    for num in s:
        if check_prime(num):
            primes.append(num)
    event= FiniteSet(*primes)
    p = probability(space, event)
    print('Sample space: {0}'.format(space))
    print('Event: {0}'.format(event))
    print('Probability of rolling a prime: {0:.5f}'.format(p))

from sympy import FiniteSet
def probability(space, event):
    return len(event)/len(space)
def check_prime(number):
    if number != 1:
        for factor in range(2, number):
            if number % factor == 0:
                return False
    else:
        return False
    return True
if __name__ == '__main__':
    space = FiniteSet(*range(1, 21))
    primes = []
    for num in space:
        if check_prime(num):
            primes.append(num)
    event= FiniteSet(*primes)
    p = probability(space, event)
    print('Sample space: {0}'.format(space))
    print('Event: {0}'.format(event))
    print('Probability of rolling a prime: {0:.5f}'.format(p))

Probability of Event A or Event B

from sympy import FiniteSet
s = FiniteSet(1, 2, 3, 4, 5, 6)
a = FiniteSet(2, 3, 5)
b = FiniteSet(1, 3, 5)
e = a.union(b)
print(len(e)/len(s))
#0.6666666666666666

Probability of Event A and Event B

from sympy import FiniteSet
s = FiniteSet(1, 2, 3, 4, 5, 6)
a = FiniteSet(2, 3, 5)
b = FiniteSet(1, 3, 5)
e = a.intersect(b)
print(len(e)/len(s))
#0.3333333333333333

Generating Random Numbers

Simulating a Die Roll

import random
print(random.randint(1, 6))

Can You Roll That Score?

   '''
   Roll a die until the total score is 20
   '''
   import matplotlib.pyplot as plt
   import random
   target_score = 20
   def roll():
       return random.randint(1, 6)
   if __name__ == '__main__':
       score = 0
       num_rolls = 0
       while score < target_score:
           die_roll = roll()
           num_rolls += 1
           print('Rolled: {0}'.format(die_roll))
           score += die_roll
       print('Score of {0} reached in {1} rolls'.format(score, num_rolls))

Is the Target Score Possible?

---

from sympy import FiniteSet
import random
def find_prob(target_score, max_rolls):
    die_sides = FiniteSet(1, 2, 3, 4, 5, 6)
    # Sample space
    s = die_sides**max_rolls
    # Find the event set
    if max_rolls > 1:
        success_rolls = []
        for elem in s:
            if sum(elem) >= target_score:
                success_rolls.append(elem)
    else:
        if target_score > 6:
            success_rolls = []
        else:
            success_rolls = []
            for roll in die_sides:
                if roll >= target_score:
                    success_rolls.append(roll)
    e = FiniteSet(*success_rolls)
    # Calculate the probability of reaching target score
    return len(e)/len(s)
if __name__ == '__main__':
    target_score = int(input('Enter the target score: '))
    max_rolls = int(input('Enter the maximum number of rolls allowed: '))
    p = find_prob(target_score, max_rolls)
    print('Probability: {0:.5f}'.format(p))

Nonuniform Random Numbers

import random
def toss():
    # 0 -> Heads, 1-> Tails
    if random.random() < 2/3:
        return 0
    else:
        return 1

'''
Simulate a fictional ATM that dispenses dollar bills
of various denominations with varying probability
'''
import random
def get_index(probability):
    c_probability = 0
    sum_probability = []
    for p in probability:
        c_probability += p
        sum_probability.append(c_probability)
    r = random.random()
    for index, sp in enumerate(sum_probability):
        if r <= sp:
            return index
    return len(probability)-1
def dispense():
    dollar_bills = [5, 10, 20, 50]
    probability = [1/6, 1/6, 1/3, 2/3]
    bill_index = get_index(probability)
    return dollar_bills[bill_index]

chapter 6: Drawing Geometric shapes and fractals

이 장에서는 fractal 그리는 연습을 해보도록 하자.

Drawing Geometric Shapes with Matplotlib’s Patches

import matplotlib.pyplot as plt
x = [1, 2, 3]
y = [1, 2, 3]
plt.plot(x, y)
#[<matplotlib.lines.Line2D object at 0x7fe822d67a20>]
plt.show()


import matplotlib.pyplot as plt
x = [1, 2, 3]
y = [1, 2, 3]
fig = plt.figure()
ax = plt.axes()
plt.plot(x, y)
#[<matplotlib.lines.Line2D object at 0x7f9bad1dcc18>]
plt.show()

Drawing a Circle

'''
Example of using matplotlib's Circle patch
'''
import matplotlib.pyplot as plt
def create_circle():
    circle = plt.Circle((0, 0), radius = 0.5)
    return circle
def show_shape(patch):
    ax = plt.gca()
    ax.add_patch(patch)
    plt.axis('scaled')
    plt.show()
if __name__ == '__main__':
    c = create_circle()
    show_shape(c)

Creating Animated Figures

'''
A growing circle
'''
from matplotlib import pyplot as plt
from matplotlib import animation
def create_circle():
    circle = plt.Circle((0, 0), 0.05)
    return circle
def update_radius(i, circle):
    circle.radius = i*0.5
    return circle,
def create_animation():
    fig = plt.gcf()
    ax = plt.axes(xlim=(-10, 10), ylim=(-10, 10))
    ax.set_aspect('equal')
    circle = create_circle()
    ax.add_patch(circle)
    anim = animation.FuncAnimation(
        fig, update_radius, fargs = (circle,), frames=30, interval=50)
    plt.title('Simple Circle Animation')
    plt.show()
if __name__ == '__main__':
    create_animation()

Animating a Projectile’s Trajectory

'''
Animate the trajectory of an object in projectile motion
'''
from matplotlib import pyplot as plt
from matplotlib import animation
import math
g = 9.8
def get_intervals(u, theta):
    flight = 2*u*math.sin(theta)/g
    intervals = []
    start = 0
    interval = 0.005
    while start < t_flight:
        intervals.append(start)
        start = start + interval
    return intervals
def update_position(i, circle, intervals, u, theta):
    t = intervals[i]
    x = u*math.cos(theta)*t
    y = u*math.sin(theta)*t - 0.5*g*t*t
    circle.center = x, y
    return circle,

def create_animation(u, theta):
    intervals = get_intervals(u, theta)
    xmin = 0
    xmax = u*math.cos(theta)*intervals[-1]
    ymin = 0
    t_max = u*math.sin(theta)/g
    ymax = u*math.sin(theta)*t_max - 0.5*g*t_max**2
    fig = plt.gcf()
    ax = plt.axes(xlim=(xmin, xmax), ylim=(ymin, ymax))
    circle = plt.Circle((xmin, ymin), 1.0)
    ax.add_patch(circle)
    anim = animation.FuncAnimation(fig, update_position,
                                    fargs=(circle, intervals, u, theta),
                                    frames=len(intervals), interval=1,
                                    repeat=False)
    plt.title('Projectile Motion')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()
if __name__ == '__main__':
    try:
        u = float(input('Enter the initial velocity (m/s): '))
        theta = float(input('Enter the angle of projection (degrees): '))
    except ValueError:
        print('You entered an invalid input')
    else:
        theta = math.radians(theta)
        create_animation(u, theta)

Drawing Fractals

Transformations of Points in a Plane

'''
Example of selecting a transformation from two equally probable
transformations
'''
import matplotlib.pyplot as plt
import random
def transformation_1(p):
    x = p[0]
    y = p[1]
    return x + 1, y - 1
def transformation_2(p):
    x = p[0]
    y = p[1]
    return x + 1, y + 1
def transform(p):
    # List of transformation functions
    transformations = [transformation_1, transformation_2]
    # Pick a random transformation function and call it
    t = random.choice(transformations)
    x, y = t(p)
    return x, y
def build_trajectory(p, n):
    x = [p[0]]
    y = [p[1]]
    for i in range(n):
        p = transform(p)
        x.append(p[0])
        y.append(p[1])
    return x, y
if __name__ == '__main__':
    # Initial point
    p = (1, 1)
    n = int(input('Enter the number of iterations: '))
    x, y = build_trajectory(p, n)
    # Plot
    plt.plot(x, y)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.show()

Drawing the Barnsley Fern

'''
Draw a Barnsley Fern
'''
import random
import matplotlib.pyplot as plt
def transformation_1(p):
    x = p[0]
    y = p[1]
    x1 = 0.85*x + 0.04*y
    y1 = -0.04*x + 0.85*y + 1.6
    return x1, y1
def transformation_2(p):
    x = p[0]
    y = p[1]
    x1 = 0.2*x - 0.26*y
    y1 = 0.23*x + 0.22*y + 1.6
    return x1, y1
def transformation_3(p):
    x = p[0]
    y = p[1]
    x1 = -0.15*x + 0.28*y
    y1 = 0.26*x + 0.24*y + 0.44
    return x1, y1
def transformation_4(p):
    x = p[0]
    y = p[1]
    x1 = 0
    y1 = 0.16*y
    return x1, y1

def get_index(probability):
    r = random.random()
    c_probability = 0
    sum_probability = []
    for p in probability:
        c_probability += p
        sum_probability.append(c_probability)
    for item, sp in enumerate(sum_probability):
        if r <= sp:
            return item
    return len(probability)-1
def transform(p):
    # List of transformation functions
    transformations = [transformation_1, transformation_2,
                       transformation_3, transformation_4]
    probability = [0.85, 0.07, 0.07, 0.01]
    # Pick a random transformation function and call it
    tindex = get_index(probability)
    t = transformations[tindex]
    x, y = t(p)
    return x, y
def draw_fern(n):
    # We start with (0, 0)
    x = [0]
    y = [0]
    x1, y1 = 0, 0
    for i in range(n):
        x1, y1 = transform((x1, y1))
        x.append(x1)
        y.append(y1)
    return x, y
if __name__ == '__main__':
    n = int(input('Enter the number of points in the Fern: '))
    x, y = draw_fern(n)
    # Plot the points
    plt.plot(x, y, 'o')
    plt.title('Fern with {0} points'.format(n))
    plt.show()

chapter 7: Sovling Calculs problems

이 장에서는 Python standard library 와 SymPy를 이용한 코딩을 해보자.

Domain and Range of a Function

An Overview of Common Mathematical Functions

import math
print(math.sin(math.pi/2))
#1.0

import sympy
sympy.sin(math.pi/2)
#1.00000000000000

from sympy import Symbol
theta = Symbol('theta')
#print(math.sin(theta) + math.sin(theta))

print(sympy.sin(theta) + sympy.sin(theta))
#2*sin(theta)


from sympy import sin, solve, Symbol
u = Symbol('u')
t = Symbol('t')
g = Symbol('g')
theta = Symbol('theta')
print(solve(u*sin(theta)-g*t, t))
#[u*sin(theta)/g]

Assumptions in SymPy

from sympy import Symbol
x = Symbol('x')
if (x+5) > 0:
    print('Do Something')
else:
    print('Do Something else')

##

x = Symbol('x', positive=True)
if (x+5) > 0:
    print('Do Something')
else:
    print('Do Something else')

Finding the Limit of Functions

from sympy import Limit, Symbol, S
x = Symbol('x')
print(Limit(1/x, x, S.Infinity))
#Limit(1/x, x, oo, dir='-')
print(l.doit())
print(Limit(1/x, x, 0, dir='-').doit())
print(Limit(1/x, x, 0, dir='+').doit())

Continuous Compound Interest

from sympy import Limit, Symbol, S
n = Symbol('n')
print(Limit((1+1/n)**n, n, S.Infinity).doit())

from sympy import Symbol, Limit, S
p = Symbol('p', positive=True)
r = Symbol('r', positive=True)
t = Symbol('t', positive=True)
print(Limit(p*(1+r/n)**(n*t), n, S.Infinity).doit())
#p*exp(r*t)

Instantaneous Rate of Change

from sympy import Symbol, Limit
t = Symbol('t')
St = 5*t**2 + 2*t + 8
t1 = Symbol('t1')
delta_t = Symbol('delta_t')
St1 = St.subs({t: t1})
St1_delta = St.subs({t: t1 + delta_t})

print(Limit((St1_delta-St1)/delta_t, delta_t, 0).doit())
#10*t1 + 2

Finding the Derivative of Functions

from sympy import Derivative, Symbol
x = Symbol('x')
f = (x**3 + x**2 + x)*(x**2+x)
print(Derivative(f, x).doit())
#(2*x + 1)*(x**3 + x**2 + x) + (x**2 + x)*(3*x**2 + 2*x + 1)

A Derivative Calculator

'''
Derivative calculator
'''
from sympy import Symbol, Derivative, sympify, pprint
from sympy.core.sympify import SympifyError
def derivative(f, var):
    var = Symbol(var)
    d = Derivative(f, var).doit()
    pprint(d)
if __name__=='__main__':
    f = input('Enter a function: ')
    var = input('Enter the variable to differentiate with respect to: ')
    try:
        f = sympify(f)
    except SympifyError:
        print('Invalid input')
    else:
        derivative(f, var)

#2*x**2 + 3*x + 1
#x

Calculating Partial Derivatives

'''
Derivative calculator
'''
from sympy import Symbol, Derivative, sympify, pprint
from sympy.core.sympify import SympifyError
def derivative(f, var):
    var = Symbol(var)
    d = Derivative(f, var).doit()
    pprint(d)
if __name__=='__main__':
    f = input('Enter a function: ')
    var = input('Enter the variable to differentiate with respect to: ')
    try:
        f = sympify(f)
    except SympifyError:
        print('Invalid input')
    else:
        derivative(f, var)

#2*x**2 + 3*x + 1
#x

#2*x*y + x*y**2
#x

Higher-Order Derivatives and Finding the Maxima and Minima

from sympy import Symbol, solve, Derivative
x = Symbol('x')
f = x**5 - 30*x**3 + 50*x
d1 = Derivative(f, x).doit()
critical_points = solve(d1)
print(critical_points)
#[-sqrt(-sqrt(71) + 9), sqrt(-sqrt(71) + 9), -sqrt(sqrt(71) + 9),
#sqrt(sqrt(71) + 9)]

A = critical_points[2]
B = critical_points[0]
C = critical_points[1]
D = critical_points[3]


d2 = Derivative(f, x, 2).doit()

print(d2.subs({x:B}).evalf())
#127.661060789073
print(d2.subs({x:C}).evalf())
#-127.661060789073
print(d2.subs({x:A}).evalf())
#-703.493179468151
print(d2.subs({x:D}).evalf())
#703.493179468151


x_min = -5
x_max = 5
print(f.subs({x:A}).evalf())
#705.959460380365
print(f.subs({x:C}).evalf())
#25.0846626340294
print(f.subs({x:x_min}).evalf())
#375.000000000000
print(f.subs({x:x_max}).evalf())
#-375.000000000000

print(f.subs({x:B}).evalf())
#-25.0846626340294
print(f.subs({x:D}).evalf())
#-705.959460380365
print(f.subs({x:x_min}).evalf())
#375.000000000000
print(f.subs({x:x_max}).evalf())
#-375.000000000000

Finding the Global Maximum Using Gradient Ascent

'''
Use gradient ascent to find the angle at which the projectile
has maximum range for a fixed velocity, 25 m/s
'''
import math
from sympy import Derivative, Symbol, sin
def grad_ascent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4
    x_old = x0
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    while abs(x_old - x_new) > epsilon:
        x_old = x_new
        x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    return x_new

def find_max_theta(R, theta):
    # Calculate the first derivative
    R1theta = Derivative(R, theta).doit()
    theta0 = 1e-3
    theta_max = grad_ascent(theta0, R1theta, theta)
    return theta_max
if __name__ == '__main__':
    g = 9.8
    # Assume initial velocity
    u = 25
    # Expression for range
    theta = Symbol('theta')
    R = u**2*sin(2*theta)/g

    theta_max = find_max_theta(R, theta)
    print('Theta: {0}'.format(math.degrees(theta_max)))
    print('Maximum Range: {0}'.format(R.subs({theta:theta_max})))

#Theta: 44.99999978475661
#Maximum Range: 63.7755102040816

A Generic Program for Gradient Ascent

'''
Use gradient ascent to find the maximum value of a
single-variable function
'''
from sympy import Derivative, Symbol, sympify
def grad_ascent(x0, f1x, x):
    epsilon = 1e-6
    step_size = 1e-4
    x_old = x0
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
while abs(x_old - x_new) > epsilon:
    x_old = x_new
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    return x_new
if __name__ == '__main__':
    f = input('Enter a function in one variable: ')
    var = input('Enter the variable to differentiate with respect to: ')
    var0 = float(input('Enter the initial value of the variable: '))
try:
    f = sympify(f)
except SympifyError:
    print('Invalid function entered')
else:
    var = Symbol(var)
    d = Derivative(f, var).doit()
    var_max = grad_ascent(var0, d, var)
    print('{0}: {1}'.format(var.name, var_max))
    print('Maximum value: {0}'.format(f.subs({var:var_max})))

#===============================================================================
# Enter a function in one variable: 25*25*sin(2*theta)/9.8
# Enter the variable to differentiate with respect to: theta
# Enter the initial value of the variable: 0.001
# theta: 0.785360029379083
# Maximum value: 63.7755100185965
#===============================================================================

#===============================================================================
# Enter a function in one variable: cos(y)
# Enter the variable to differentiate with respect to: y
# Enter the initial value of the variable: 0.01
# y: 0.00999900001666658
# Maximum value: 0.999950010415832
#===============================================================================

#===============================================================================
# Enter a function in one variable: cos(y) + k
# Enter the variable to differentiate with respect to: y
# Enter the initial value of the variable: 0.01
# y: 0.00999900001666658
# Maximum value: k + 0.999950010415832
#===============================================================================

#===============================================================================
# Enter a function in one variable: x**5 - 30*x**3 + 50*x
# Enter the variable to differentiate with respect to: x
# Enter the initial value of the variable: -2
# x: -4.17445116397103
# Maximum value: 705.959460322318
#===============================================================================

#===============================================================================
# Enter a function in one variable: x**5 - 30*x**3 + 50*x
# Enter the variable to differentiate with respect to: x
# Enter the initial value of the variable: 0.5
# x: 0.757452532565767
# Maximum value: 25.0846622605419
#===============================================================================

The Role of the Step Size and Epsilon

'''
Use gradient ascent to find the maximum value of a
single-variable function. This also checks for the existence
of a solution for the equation f'(x)=0.
'''
from sympy import Derivative, Symbol, sympify, solve
def grad_ascent(x0, f1x, x):
    # Check if f1x=0 has a solution
    if not solve(f1x):
        print('Cannot continue, solution for {0}=0 does not exist'.format(f1x))
        return
    epsilon = 1e-6
    step_size = 1e-4
    x_old = x0
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
while abs(x_old - x_new) > epsilon:
    x_old = x_new
    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()
    return x_new

if __name__ == '__main__':
    f = input('Enter a function in one variable: ')
    var = input('Enter the variable to differentiate with respect to: ')
    var0 = float(input('Enter the initial value of the variable: '))
    try:
        f = sympify(f)
    except SympifyError:
        print('Invalid function entered')
    else:
        var = Symbol(var)
        d = Derivative(f, var).doit()
        var_max = grad_ascent(var0, d, var)
    if var_max:
        print('{0}: {1}'.format(var.name, var_max))
        print('Maximum value: {0}'.format(f.subs({var:var_max})))

#===============================================================================
# Enter a function in one variable: log(x)
# Enter the variable to differentiate with respect to: x
# Enter the initial value of the variable: 0.1
# Cannot continue, solution for 1/x=0 does not exist
#===============================================================================

Finding the Integrals of Functions

from sympy import Integral, Symbol



if __name__=="__main__":
    x = Symbol('x')
    k = Symbol('k')
    print(Integral(k*x, x))
    #Integral(k*x, x)

    print(Integral(k*x, x).doit())
    #k*x**2/2
    print(Integral(k*x, (x, 0, 2)).doit())
    #2*k


from sympy import Integral, Symbol
x = Symbol('x')
print(Integral(x, (x, 2, 4)).doit())

Probability Density Functions

from sympy import Symbol, exp, sqrt, pi, Integral
x = Symbol('x')
p = exp(-(x - 10)**2/2)/sqrt(2*pi)
print(Integral(p, (x, 11, 12)).doit().evalf())
#0.135905121983278


from sympy import Symbol, exp, sqrt, pi, Integral, S
x = Symbol('x')
p = exp(-(x – 10)**2/2)/sqrt(2*pi)
print(Integral(p, (x, S.NegativeInfinity, S.Infinity)).doit().evalf())
#1.00000000000000