Developers Heaven Forum
Desktop Programming => HLA => Topic started by: David on January 14, 2013, 12:56:16 AM
-
Below, please find the Python Demo Programs linked to here, from the new thread ...
Beginning Calculus
To go to Beginning Calculus, click here ...
http://developers-heaven.net/forum/index.php?topic=2601.0
To go to a Beginner's Level introduction to the Python version 3 programming language ...
http://developers-heaven.net/forum/index.php?topic=46.0
This first demo is an 'updated version' of Archimedes method to calculate the value for the constant pi
Select/Copy the following program, coded in the Python computer language, and paste that code into a Text Editor like Notepad ... then 'SaveAs' that file with file name my_piArchimedes.py
Making sure, that you have a 3.x version of the freely available Python Interpreter installed on your PC, then, you can 'click on' ... i.e. 'run' that .py file you just saved.
Note: Python program file names end with .py
# my_piArchimedes.py # this ver. 2013-01-20 #
# This Python program needs Python version 3 (or newer version)
# installed on your PC ... to run ok
header1 = \
'''A demo program in 'Beginning Calculus' by dwzavitz
Find perimeter/2 of inscribed regular polygon for circle radius = 1 unit
using here an ... 'Updated Archimedes Method' ...
that cal's sin(t/2) from sin(t) as the number of sides double each loop
Uses sin(t/2) = sqrt( (1 - sqrt(1 - (sin(t))^2) ) / 2 ) ... (see below)
for regular inscribed polygon ... (for 'UNIT CIRCLE')
_
|_|_r=1, n=4, t/2 = 45 deg's, sin(t/2) = 1/sqrt(2), pi = 4/sqrt(2) = 2.824
_
/ \_r=1, n=6, t/2 = 30 deg's, sin(t/2) = 0.5, s = 1, pi = n*(0.5) = 3
\_/
'''
header2 = \
'''Recall:
e^(i*2*t) = cos(2*t) + i*sin(2*t)
e^(i*2*t) = e^(i*t)^2 = (cos(t) + i*sin(t))^2
= (cos(t))^2 - (sin(t))^2) + i*2*sin(t)cos(t)
And taking 'non-i' parts to be equal ...
cos(2*t) = (cos(t))^2 - (sin(t))^2)
So cos(2*t) = 1 - 2*(sin(t))^2 ... Recall (sin(x))^2 + (cos(x))^2 = 1
So sin(t) = sqrt( (1 - cos(2*t) ) / 2 )
And sin(t/2) = sqrt( (1 - cos(t) ) / 2 )
But since ... ... cos(t) = sqrt(1 - (sin(t))^2)
Now sin(t/2) = sqrt( (1 - sqrt(1 - (sin(t))^2) ) / 2 )
'''
from math import * # make pi (and all math library) available in global namespace #
# above '*' char indicates to import all math functions, etc...
# and so now ... don't need 'math.xxx' math namespace name at front of 'xxxdef Perimeter( n, sin_th ):
def halfPerimeter( n, sin_th ): ## Note: (n/2 - 2*sin_th) = n * sin_th ##
return n * sin_th
print( header1 )
input( "Press 'Enter' to continue ... " )
print( header2 )
input( "Press 'Enter' to continue ... " )
print( "Loop# Sides sin(th/2) " + \
"Half Perimeter -> Area ('Unit Circle')" )
#n = 4 ## n is number of sides to start (of inscribed regular polygon) ##
n = 6 ## Note: 6 sides was Archimedes starting number ##
#sin_th = 2**0.5/2 # start value: 4 sides and t = 90, t/2 = 45 deg's
sin_th = 0.5 # start value: 6 sides and t = 60, t/2 = 30 ...
P = 0
count = 0
while True: # loop forever until break condition below is realized #
count += 1
P_old = P
P = halfPerimeter( n, sin_th )
if count % 5 == 1:
print( '{:2d} {:11,d} {:.15e} {:.16f}'.format(count,n,sin_th, P) )
if P - P_old < 10**(-16): # i.e. when reach limit of 'float precision'
break
# double the number of sides to get ready for next loop
n = 2*n
# now update the sine of the next 1/2 angle when sides doubled
#sin(t/2) = sqrt( (1 - sqrt(1 - (sin(t))^2) ) / 2 )
#sin_th = sqrt( (1 - sqrt(1-sin_th**2) / 2 ) #updating to sin(th/2) here#
## MUCH better sig dig result in PC if use this form of above ##
## by multiplying top and bottom of fraction by sqrt( (1+sqrt(1-sin_th**2)) ) ##
sin_th = sin_th / 2**0.5 /sqrt( 1 + sqrt( 1- sin_th**2 ) )
print() # print a blank line ...
print( '%.16f' % pi, '<- pi (the actual value, rounded to 16 decimal places)' )
print( '%.16f %.16f <- P_old' % (P_old, P_old/pi) )
print( '%.16f %.16f <- P' % (P, P/pi) )
print( '%.16f <- (P - P_old) compared to 10**-16 = %.16f' % (P - P_old, 10**-16) )
input( "\nPress 'Enter' to continue/exit ... " )
# program output ... when run ...
'''
Press 'Enter' to continue ...
Loop# Sides sin(th/2) Half Perimeter -> Area ('Unit Circle')
1 6 5.000000000000000e-01 3.0000000000000000
6 192 1.636173162648678e-02 3.1414524722854611
11 6,144 5.113269070136972e-04 3.1415925166921559
16 196,608 1.597896653979543e-05 3.1415926534561005
21 6,291,456 4.993427043898358e-07 3.1415926535896590
26 201,326,592 1.560445951218302e-08 3.1415926535897896
3.1415926535897931 <- pi (the actual value, rounded to 16 decimal places)
3.1415926535897896 0.9999999999999989 <- P_old
3.1415926535897896 0.9999999999999989 <- P
0.0000000000000000 <- (P - P_old) compared to 10**-16 = 0.0000000000000001
Press 'Enter' to continue/exit ...
'''
This next demo finds perimeter/2 of BOTH inscribed AND exscribed regular polygons ... for a circle with radius = 1 unit
(using here an ... 'Updated Archimedes ***SQUEEZE*** Method')
# my_piArchimedes2.py # this ver. 2013-01-20 #
# This Python program needs Python version 3 (or newer version)
# installed on your PC ... to run ok
header1 = \
'''A demo program in 'Beginning Calculus' by dwzavitz
Find perimeter/2 of inscribed regular polygon for circle radius = 1 unit
AND perimeter/2 of EXscribed regular polygon ...
using here an ... 'Updated Archimedes ***SQUEEZE*** Method' ...
that cal's sin(t/2) from sin(t) as the number of sides double each loop
Uses sin(t/2) = sin(t)/2^0.5/sqrt( 1+sqrt(1-(sin(t))^2) ) ... (see below)
And tan(t/2) = tan(t) / (sqrt(1+(tan(t))^2)) + 1) ... (see below)
for regular inscribed polygon ... (for 'UNIT CIRCLE')
_
|_|_r=1, n=4, t/2 = 45 deg's, sin(t/2) = 1/sqrt(2), pi = 4/sqrt(2) = 2.824...
_
/ \_r=1, n=6, t/2 = 30 deg's, sin(t/2) = 0.5, s = 1, pi = n*(0.5) = 3
\_/ for EXscribed hexagon tan(t/2) = 1/sqrt(3), pi = n/sqrt(3) = 3.464...
'''
header2 = \
'''Recall:
e^(i*2*t) = cos(2*t) + i*sin(2*t)
e^(i*2*t) = e^(i*t)^2 = (cos(t) + i*sin(t))^2
= (cos(t))^2 - (sin(t))^2) + i*2*sin(t)cos(t)
And taking 'non-i' parts to be equal ...
cos(2*t) = (cos(t))^2 - (sin(t))^2)
So cos(2*t) = 1 - 2*(sin(t))^2 ... Recall (sin(x))^2 + (cos(x))^2 = 1
So sin(t) = sqrt( (1 - cos(2*t) ) / 2 )
And sin(t/2) = sqrt( (1 - cos(t) ) / 2 )
But since ... ... cos(t) = sqrt(1 - (sin(t))^2)
(So for inscribed regular polygon ...)
sin(t/2) = sqrt( (1 - sqrt(1 - (sin(t))^2) ) / 2 )
(And for excribed regular polygon ...)
tan(t) = 2*tan(t/2) / ( 1 - (tan(t/2))^2 ) ... implies
c = 2*x / (1-x^2) ...substituting in c = tan(t), x = tan(t/2) implies
x = ( sqrt(1+c^c) - 1 ) / c
= c / ( sqrt(1+c^2) +1 ) ... implies ...
tan(t/2) = tan(t) / (sqrt(1+(tan(t))^2)) + 1)
'''
# make pi, sqrt (and all math library) available in global namespace #
from math import *
# above '*' char indicates to import all math functions, etc...
# and so now ... don't need 'math.xxx' math namespace name at front of 'xxxdef Perimeter( n, sin_th ):
def halfPerimeter( n, sin_th ): ## Note: (n/2 - 2*sin_th) = n * sin_th ##
return n * sin_th
print( header1 )
input( "Press 'Enter' to continue ... " )
print( header2 )
input( "Press 'Enter' to continue ... " )
print()
print( "Loop# Sides sin(th/2) / tan(th/2) " + \
"Half Perimeter -> Area ('Unit Circle')" )
#n = 4 ## n is number of sides to start (of inscribed regular polygon) ##
n = 6 ## Note: 6 sides was Archimedes starting number ##
#sin_th = 2**0.5/2 # start value: 4 sides and t = 90, t/2 = 45 deg's
sin_th = 0.5 # start value: 6 sides and t = 60, t/2 = 30 ...
#tan_th = 1 # for starting with exscribed square #
tan_th = 3**(-0.5) # for starting exscribed hexagon #
P = 0
count = 0
while True: # loop forever until break condition below is realized #
count += 1
P_old = P
P = halfPerimeter( n, sin_th )
Ptan = tan_th * n
if count % 5 == 1:
print( '{:2d} {:11,d} {:.15e} {:.16f}'.format(count,n,sin_th, P) )
print( '{:2d} {:11,d} {:.15e} {:.16f}'.format(count,n,tan_th, Ptan) )
print()
if P - P_old < 10**(-16): # i.e. when reach limit of 'float precision'
break
# double the number of sides to get ready for next loop
n = 2*n
# now update the sine of the next 1/2 angle when sides doubled
#sin(t/2) = sqrt( (1 - sqrt(1 - (sin(t))^2) ) / 2 )
#sin_th = sqrt( (1 - sqrt(1-sin_th**2) / 2 ) #updating to sin(th/2) here#
## MUCH better sig dig result in PC if use this form of above ##
## by multiplying top and bottom of fraction by sqrt( (1+sqrt(1-sin_th**2)) ) ##
sin_th = sin_th / 2**0.5 /sqrt( 1 + sqrt( 1- sin_th**2 ) )
tan_th = tan_th / (sqrt(1+tan_th**2) + 1)
#print() # print a blank line ...
print( '{:.16f} <- average of inscribed and exscribed ...'.format( (P+Ptan)/2 ) )
print( '%.16f' % pi, '<- pi (the actual value, rounded to 16 decimal places)' )
input( "\nPress 'Enter' to continue/exit ... " )
# program output ... when run ...
'''
Press 'Enter' to continue ...
Loop# Sides sin(th/2) / tan(th/2) Half Perimeter -> Area ('Unit Circle')
1 6 5.000000000000000e-01 3.0000000000000000
1 6 5.773502691896257e-01 3.4641016151377544
6 192 1.636173162648678e-02 3.1414524722854611
6 192 1.636392213531158e-02 3.1418730499798242
11 6,144 5.113269070136972e-04 3.1415925166921559
11 6,144 5.113269738582516e-04 3.1415929273850978
16 196,608 1.597896653979543e-05 3.1415926534561005
16 196,608 1.597896654183539e-05 3.1415926538571730
21 6,291,456 4.993427043898358e-07 3.1415926535896590
21 6,291,456 4.993427043898989e-07 3.1415926535900560
26 201,326,592 1.560445951218302e-08 3.1415926535897896
26 201,326,592 1.560445951218305e-08 3.1415926535897962
3.1415926535897931 <- average of inscribed and exscribed ...
3.1415926535897931 <- pi (the actual value, rounded to 16 decimal places)
Press 'Enter' to continue/exit ...
'''
This next Python code, takes good advantage of Pythons built in unlimited-long integers, to find pi to up 100,000 decimal places ...
Enjoy ...
# my_pi100_000places.py # this ver. 2013-01-11 #
# This file to be run under Python ver. 3 #
header = \
'''A demo program in 'Beginning Calculus' by dwzavitz
to find pi up to 100,000 decimal places
'''
# acotx = 1/x - 1/3x**3 + 1/5x**5 - 1/7x**7 + ... #
def acot( x, unity ):
ssum = xpower = unity // x
n = 3
sign = -1
count = 0
while 1:
xpower = xpower // (x*x)
term = xpower // n
if not term:
break
ssum += sign*term
sign = -sign
n += 2
count += 1
'''
if count == 100000:
print( "count = ", count )
break
'''
print( "count = {:,d}".format( count ) )
#input( "\nPress 'Enter' to continue ... " )
return ssum
def mypi( digits ):
unity = 10**( digits + 10 )
# PI/4 <= 4*ACOT5 - ACOT239 #
return 4 * ( 4 * acot( 5, unity ) - acot(239, unity) ) // 10**10
print( header )
done = False
while not done:
s = 'Num of xxx decimal places to find in pi = 3.xxx (20 to 100000): '
try:
dig = int(input( s ))
done = True
if dig > 100000:
done = False
raise
elif dig >= 30000:
done = False
print('{:,d}'.format(dig),"could take several seconds ... ",end = '')
if input( "Ok (y/n) ? " ) == 'y' :
done = True
if dig < 20:
dig = 20
print( "Doing min of 20 decimal places ..." )
except:
print( "Invalid ... integer input <= 100,000 only please ... " )
import datetime
t1 = datetime.datetime.now()
tostr = str(mypi(dig)) # call mypi and convert num to str ...
print( "len(tostr)", '{:,d}'.format( len(tostr) ) )
import sys
outSave = sys.stdout
fout = open( "PI.txt", "w" )
sys.stdout = fout
print( "mypi(", dig, ")=\n", tostr[:1] + '.', sep = '' )
count = 0
for c in tostr[1:]:
count += 1
print( c, sep='', end='' )
if count == 100:
print()
count = 0
#print()
#print( '3.\n' + '1234567890'*(dig//10) )
sys.stdout = outSave
fout.close()
t2 = datetime.datetime.now()
print( "Processing time was (time2 - time1) =", t2-t1 )
import os
os.system( 'notepad ' + 'PI.txt' )
#input( "\nPress 'Enter' to continue/exit ... " )
-
This next little demo, coded in Python 3, shows a way to find the AREA under a curve ...
The integral here summed is ...
the area under the curve of 1/t
for t = 1 to t = x
( i.e. x >= e ... because the program stops when the AREA under the curve becomes >= 1 )
# ln_x_by_area.py * this ver 2013-01-17 #
header = \
'''A demo program in 'Beginning Calculus' by dwzavitz
Finds e by finding area under graph of 1/t
as t goes from 1 to x
stops at t = e ... when area = sum_of_slices >= 1.0
In a loop, repeats above,
with slices 1/10 width of former width
stop this (outer) loop when number of slices n > 10**8
'''
print( header )
input( "Press 'Enter' to continue ... " )
print()
import sys
def mySum( n ):
area = 0.0
w = 3 - 1 # max width of interval
d = w/n # increment size
t = 1-d # get ready to increment
print( '{:,d}'.format(n), '(max) slices sum to ', end = '' )
sys.stdout.flush()
for looop in range( n ):
t += d
y = 1/t
area += y*d
if area >= 1.0:
break
return area, t
n = 100 # start here with length 2 divided into 100 equal parts
count = 0
while True:
area, t = mySum( n )
s1 = '{:.' + str(len(str(n))-1) + 'f}'
s2 = '{:.' + str(len(str(n))-2) + 'f}'
count += 1
print( s1.format( area ), 'at estimate({}, e) = t ='.format(count), \
s2.format(t) )
n *= 10
if n == 10**7 :
print( "\nDoing 2nd last loop now ... this will take a few sec's" )
if n == 10**8 :
print( "\nDoing last loop now ... this will take several sec's" )
if n > 10**8:
break
input( "\nPress 'Enter' to continue/exit ... " )
# output #
'''
A demo program in 'Beginning Calculus' by dwzavitz
Finds e by finding area under graph of 1/t
as t goes from 1 to x
stops at t = e ... when area = sum_of_slices >= 1.0
in a loop, repeats above,
with slices 1/10 width of former width
stop this (outer) loop when number of slices n > 10**8
100 slices sum to 1.01 at estimate(1, e) = t = 2.7
1,000 slices sum to 1.001 at estimate(2, e) = t = 2.72
10,000 slices sum to 1.0000 at estimate(3, e) = t = 2.718
100,000 slices sum to 1.00001 at estimate(4, e) = t = 2.7183
1,000,000 slices sum to 1.000001 at estimate(5, e) = t = 2.71828
Doing 2nd last loop now ... this will take a few sec's
10,000,000 slices sum to 1.0000001 at estimate(6, e) = t = 2.718282
Doing last loop now ... this will take several sec's
100,000,000 slices sum to 1.00000000 at estimate(7, e) = t = 2.7182818
Press 'Enter' to continue/exit ...
'''
A 'fast' version of the above that utilizes the property ...
Area under curve of 1/t from 1 to a^n = n * (Area under curve of 1/t from 1 to a)
i.e. ln(a^n) = n * ln(a)
# ln_x_by_areaFAST.py * this ver 2013-01-17 #
header = \
'''A demo program in 'Beginning Calculus' by dwzavitz
Finds e by finding area under graph of 1/t
as t goes from 1 to x
stops at t = e ... when area = sum_of_slices >= 1.0
in a loop, repeats above,
with slices 1/10 width of former width
stop this (outer) loop when number of slices n > 10**8
'''
print( header )
input( "Press 'Enter' to continue ... " )
print()
import sys
import math
log = math.log10
def mySum( n, stop ):
area = 0.0
w = 3 - 1 # max width of interval
d = w/n # increment size
print( '{:,d}'.format(n), '(max) slices sum to ', end = '' )
sys.stdout.flush()
t = 1 + d
mult = t
y = 1/t
area = d # but use t = 1 # area = y*d
loops = 1
'''
while True:
loops += 1
t *= mult
if t >= 2.5:
break
'''
# if mult^n >= 2.5
# then n*log10(mult) = log10(2.5)
# and n = int(math.ceil(...))
loops = int(math.ceil(log(stop)/log(mult)))
'''
print( 'int(math.ceil(log(2.5)/log(mult))) =', \
int(math.ceil(loops)), 'loops =', loops )
'''
area *= loops
t = t**loops
loops2 = 0
while area < 1:
loops2 += 1
t += d
y = 1/t
area += y*d
#print( 'loops2 =', loops2 )
return area, t, loops, loops2
n = 100 # start here with length 2 divided into 100 equal parts
count = 0
s = '{:,d}'
s3 = 'Actual slices = ' + s + ' + ' + s + ' = ' + s
stop = 2.5
while True:
area, t, loops, loops2 = mySum( n, stop )
s1 = '{:.' + str(len(str(n))-1) + 'f}'
s2 = '{:.' + str(len(str(n))-2) + 'f}'
count += 1
print( s1.format( area ), 'at estimate({}, e) = t ='.format(count), \
s2.format(t) )
print( s3.format(loops, loops2, loops+loops2 ) )
n *= 10
if n > 10**8:
break
stop = t
input( "\nPress 'Enter' to continue/exit ... " )
# output #
'''
A demo program in 'Beginning Calculus' by dwzavitz
Finds e by finding area under graph of 1/t
as t goes from 1 to x
stops at t = e ... when area = sum_of_slices >= 1.0
in a loop, repeats above,
with slices 1/10 width of former width
stop this (outer) loop when number of slices n > 10**8
Press 'Enter' to continue ...
100 (max) slices sum to 1.00 at estimate(1, e) = t = 2.7
Actual slices = 47 + 8 = 55
1,000 (max) slices sum to 1.001 at estimate(2, e) = t = 2.72
Actual slices = 497 + 9 = 506
10,000 (max) slices sum to 1.0000 at estimate(3, e) = t = 2.718
Actual slices = 4,999 + 3 = 5,002
100,000 (max) slices sum to 1.00001 at estimate(4, e) = t = 2.7183
Actual slices = 49,997 + 9 = 50,006
1,000,000 (max) slices sum to 1.000000 at estimate(5, e) = t = 2.71828
Actual slices = 499,999 + 3 = 500,002
10,000,000 (max) slices sum to 1.0000001 at estimate(6, e) = t = 2.718282
Actual slices = 4,999,997 + 9 = 5,000,006
100,000,000 (max) slices sum to 1.00000000 at estimate(7, e) = t = 2.7182818
Actual slices = 49,999,999 + 3 = 50,000,002
Press 'Enter' to continue/exit ...
'''
-
Find e^x up to 200,000 places of decimals ... (makes good use of Pythons unlimited-long integers) ...
First ... just e
Then following ... e^x
# my_e200_000places.py # this ver. 2013-01-13 #
# This file to be run under Python ver. 3 #
header = \
'''A demo program in 'Beginning Calculus' by dwzavitz
to find e up to 200,000 decimal places
'''
def getE(x, shift):
term = 10**(shift+10)
ex = divisor = 0
while True:
exOld = ex
ex += term
divisor += 1
#term = term * x // divisor
term = term // divisor
if exOld == ex:
break;
return ex // 10**10 # trim back off 10 extra places shifted left
def takeIn():
done = False
n = 1
while not done:
try:
#n = int(input("Input 'n' to find e**n : "))
d = int(input("Decimal digits desired (1,000 to 200,000): "))
if d < 1000:
d = 1000
print( "(Doing min of 1,000)", end = ' ' )
if d >= 50000:
print( 'This will take a few seconds ...' , end = ' ' )
if d > 200000:
print( "(Doing max of 200,000)", end = ' ' )
d = 200000
done = True
except:
print( 'Invalid entry ... numbers only please ...' )
return n, d
print( header )
import sys
import os
while 1:
n, d = takeIn()
print()
outSave = sys.stdout
sys.stdout = open( "my_e.txt", "w" )
s = str( getE(n, d) )
print( len(s), 'digits in total ... ' )
print( s[0] + '.' )
s = s[1:]
for n, c in enumerate(s):
print( c, end = '' )
if (n+1) % 100 == 0: print()
sys.stdout = outSave
os.system( 'notepad ' + 'my_e.txt' )
reply = input( "More (y/n) " )
if( reply and reply in "Nono" ):
break
Now ... e^x ...
# my_eX200_000places.py # this ver. 2013-01-13 #
# This file to be run under Python ver. 3 #
header = \
'''A demo program in 'Beginning Calculus' by dwzavitz
to find e^x up to 200,000 decimal places
'''
import math
def getE(x, shift):
places = math.ceil(x * math.log10(math.e))
if places > 0 :
term = 10**(shift+10+places)
else:
term = 10**(shift+10+17)
ex = divisor = 0
toInt = 10**17
n = int( x * toInt )
while True:
exOld = ex
ex += term
divisor += 1
term = term * n // divisor
term = term // toInt
if exOld == ex:
break;
s = str(ex)
sLen = len(s)
if places > 0:
#return ex // 10**10 # trim back off (10+places) \
# extra places shifted left
return s[0:-(10+places)]
#return ex // 10**17 # trim back off 17 extra places shifted left
return s[0:-17]
def takeIn():
done = False
while not done:
try:
x = eval(input("Input 'x' to find e**x: "))
if x < 10**-15 or x > 700 :
print( "10^-15 <= x <= 700 is the "
"allowed range here.", end = ' ' )
raise
d = int(input("Decimal digits desired (1,000 to 200,000): "))
if d < 1000:
print( "(Doing min of 1,000)", end = ' ' )
d = 1000
if d >= 50000:
print( 'This will take a few seconds ...' , end = ' ' )
if d > 200000:
print( "(Doing max of 200,000)", end = ' ' )
d = 200000
done = True
except:
print( 'Error!!! Try again ...' )
return x, d
print( header )
import sys
import os
while 1:
x, d = takeIn()
print()
outSave = sys.stdout
sys.stdout = open( "my_eX.txt", "w" )
#s = str( getE(x, d))
s = getE( x, d )
lenS = len(s)
decPlaceIn = lenS - d
print( math.exp(x), '= math.exp({}) so can ' \
'compare digits and exponent.'.format(x) )
print( lenS, 'digits in total ... ' )
#print( s[0:decPlaceIn] + '.' )
for n, c in enumerate( s[0:decPlaceIn] ):
print( c, end = '' )
if (n+1) % 100 == 0: print()
print( '.' )
# print( s[0] + '.' )
s = s[decPlaceIn:]
for n, c in enumerate(s):
print( c, end = '' )
if (n+1) % 100 == 0: print()
sys.stdout = outSave
os.system( 'notepad ' + 'my_eX.txt' )
reply = input( "More (y/n) " )
if( reply and reply in "Nono" ):
break
-
A little demo to show that, as x increases, then e^x exceeds x^n for all finite n ...
# table_eX_vs_xN.py # this ver. 2013-01-22 #
import math
exp = math.exp
def signLog( x ):
if x == 0:
return '-'
elif math.log10(x) >= 0:
return '+'
else:
return '-'
def pause( msg ):
input( "Press 'Enter' to continue{} ... ".format( msg) )
for n in range(2, 12):
print( " x \t x^{} \t\t e^x \t\t x^{}/e^x".format( n, n ), \
' \t signLog(x^{}/e^x)'.format(n) )
for x in range(0, 10):
xx = x/10
b = exp(xx)
a = xx**n
c = a/b
print( '{:>2} {:>16e} {:>16e} {:>16e}'. format( xx, a, b, c), \
' ', signLog(c) )
pause( '' )
for x in range(1, 60):
b = exp(x)
a = x**n
c = a/b
print( '{:>2} {:>16e} {:>16e} {:>16e}'. format( x, a, b, c), \
' ', signLog(c) )
if x % 20 == 0:
pause( '' )
print()
pause( '/exit')
A little demo to show that, as x increases, then 2^x exceeds x^n for all finite n ...
# table_bX_vs_xN.py # this ver. 2013-01-20 #
import math
##exp = math.exp
base = 2
def signLog( x ):
if x == 0:
return '-'
elif math.log10(x) >= 0:
return '+'
else:
return '-'
def pause( msg ):
input( "Press 'Enter' to continue{} ... ".format( msg) )
for n in range(2, 12):
print( " x \t x^{} \t\t {}^x \t\t x^{}/{}^x".format(n, base, n, base), \
' \t signLog(x^{}/{}^x)'.format(n, base) )
for x in range(0, 10):
xx = x/10
#b = exp(xx)
b = base**xx
a = xx**n
c = a/b
print( '{:>2} {:>16e} {:>16e} {:>16e}'. format( xx, a, b, c), \
' ', signLog(c) )
pause( '' )
for x in range(1, 70):
#b = exp(x)
b = base**x
a = x**n
c = a/b
print( '{:>2} {:>16e} {:>16e} {:>16e}'. format( x, a, b, c), \
' ', signLog(c) )
if x % 20 == 0:
pause( '' )
print()
pause( '/exit')
-
Did you ever want to find the roots of some polynomial ?
Here is 'Newton's Method' ... to the rescue.
Roots by Newton's (iterative) method of successive (improved) approximations
Xn = Xp - f(Xp)/f'(Xp) i.e. Xnext = Xprev - f'(Xprev)/f(Xprev)
.
./ |f(Xp)
./ |
./ |
. /__|
. Xn Xp Method to find Xn is derived from:
f(Xp)/(Xp-Xn) = rise/run = slope = f'(Xp)
# newtonsRoots.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
testData = [-10, -1.1, -0.1, 0.1, 1.1, 10]
'''
testData = [-10, -1.1, -0.1, 0.1, 1.1, 10] # x[i] to use in testData
print( header )
print( '\nloops \t\t prev x \t new x')
def f( x ):
return x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
def df( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if 1 - err <= (a/b) <= 1 + err:
return True
return False
def inList( v, lst ):
for item in lst:
if isEqual( v, item ): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
s = '{:<3} {:<20} {:<20}'
print( s.format(count, x0, x1) )
print( "Check if f({}) = {} = 0 (or nearly so?)".format( x1, f(x1) ) )
if not inList( x1, roots ):
roots.append( x1 )
input("Press 'Enter' to continue ... ")
print()
print("\n\nNow ... finding roots of df = f2 (to find max and min) ...")
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual( v, item[1] ): return True
return False
def f2( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def df2( x ):
return 20*x**3 - 48*x**2 - 30*x + 40
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1): #math.fabs(x1-x0) > 10**-15: # x0/x1 != 1:
count += 1
x0 = x1
x1 = getx2(x1)
s = '{:<3} {:<20} {:<20}'
print( s.format(count, x0, x1) )
print( "Check if f2({}) = {} = 0 (or nearly so?)".format( x1, f2(x1) ) )
tmp = df2(x1)
s2 = '{:.2f}'
print( "Value of df2({}) = {}".format(x1, tmp), \
maxMin(tmp), s2.format(x1) )
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
numPoints = 50
intervalWidth = 9
end = intervalWidth/2
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints, s.format(-end, end-dw) )
points = []
# create table #
for p in range( numPoints ):
v = -end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print( roots )
print( "MaxMin:" )
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
The above, a little more polished ...
# newtonsRoots5.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
testData = [x/10 for x in range(-101, 102)]
'''
testData = [x/10 for x in range(-101, 102)] # x[i] to use in testData #
print( header )
def f( x ):
return x**5 - 4*x**4 - 5*x**3 + 20*x**2 + 4*x - 16
def df( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if 1 - err <= (a/b) <= 1 + err:
return True
return False
def inList( v, lst ):
for item in lst:
if isEqual( v, item ): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
if not inList( x1, roots ):
roots.append( x1 )
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual( v, item[1] ): return True
return False
def f2( x ):
return 5*x**4 - 16*x**3 - 15*x**2 + 40*x + 4
def df2( x ):
return 20*x**3 - 48*x**2 - 30*x + 40
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx2(x1)
tmp = df2(x1)
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
roots.sort()
numPoints = 50
intervalWidth = roots[-1]+1 - roots[0]
end = roots[0] - 1
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints+5, s.format(end, end+(numPoints+4)*dw) )
points = []
# create table #
for p in range( numPoints+5 ):
v = end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print('[', end = '' )
for i, r in enumerate(roots):
print( '{:.3f}'.format(r), end = '' )
if i < len(roots)-1: print(', ', end = '' )
else: print(']')
print( "\nMaxMin:" )
from operator import itemgetter
extremes.sort(key=itemgetter(1))
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
The above, but a more general solution ...
# newtonsRoots6.py # this ver. 2013-02-01 #
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
Finds roots and max's and min's of polynomial using Newton's method
Example used here is:
f(x) = x**6 - 4*x**5 - 5*x**4 + 20*x**3 + 4*x**2 - 16*x
testData = [x/10 for x in range(-1001, 1002)]
'''
testData = [x/10 for x in range(-1001, 1002)] # x[i] to use in testData #
print( header )
def f( x ):
return x**6 - 4*x**5 - 5*x**4 + 20*x**3 + 4*x**2 - 16*x
def df( x ):
return 6*x**5 - 20*x**4 - 20*x**3 + 60*x**2 + 8*x -16
def getx1( x ):
return x - f(x)/df(x)
import math
err = 10**-15
def isEqual( a, b ):
if b != 0:
if 1 - err <= (a/b) <= 1 + err:
return True
return False
else:
if abs(b-a) < err: return True
else: return False
def inList( v, lst ):
for item in lst:
if isEqual(v, item): return True
return False
roots = []
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx1(x1)
if not inList( x1, roots ):
roots.append( x1 )
extremes = []
def inList( v, lst ):
for item in lst:
if isEqual(v, item[1]): return True
return False
def f2( x ):
return 6*x**5 - 20*x**4 - 20*x**3 + 60*x**2 + 8*x -16
def df2( x ):
return 30*x**4 - 80*x**3 - 60*x**2 + 120*x + 8
def getx2( x ):
return x - f2(x)/df2(x)
def maxMin( x ):
if x < 0: return "=> a Max at"
elif x > 0: return "=> a Min at"
else: return "=> an inflection at"
for x1 in testData:
if x1 == 0: continue
x0 = 0
count = 0
while not isEqual(x0,x1):
count += 1
x0 = x1
x1 = getx2(x1)
tmp = df2(x1)
if not inList( x1, extremes ):
extremes.append( [maxMin(tmp), x1, f(x1)] )
input("Press 'Enter' to continue ... ")
print()
roots.sort()
numPoints = 50
intervalWidth = roots[-1]+1 - roots[0]
end = roots[0] - 1
dw = intervalWidth/numPoints
s = "points (x,y), on curve y = f(x), with x in [{:.2f},{:.2f}] ...\n"
print( "Table of", numPoints+5, s.format(end, end+(numPoints+4)*dw) )
points = []
# create table #
for p in range( numPoints+5 ):
v = end + p*dw
points.append( [v, f(v)] )
def getSign(x):
if x<0: return 0
return 1
# show table #
pOld = points[0][1]
n = 0
for p in points:
print( '({:5.2f},{:7.1f})'.format(p[0],p[1]), end = '' )
n += 1
if n == 5:
print()
n = 0
else:
print( ' ', end = '' )
if getSign(p[1]) != getSign(pOld): print( "<- sign changed\n" )
pOld = p[1]
print( "\nRoots:" )
print('[', end = '' )
for i, r in enumerate(roots):
print( '{:.3f}'.format(r), end = '' )
if i < len(roots)-1: print(', ', end = '' )
else: print(']')
print( "\nMaxMin:" )
from operator import itemgetter
extremes.sort(key=itemgetter(1))
for mm in extremes:
print( mm )
input("Press 'Enter' to continue ... ")
Just find nth root of some value (using Newton's Method) ...
# newtonsRootsNth.py # this ver. 2013-02-01 #
# nth root of value demo #
err = 10**-15
header = \
'''Demo program linked from 'Beginning Calculus' by dwzavitz
1. Input positive value to find root ...
2. Input integer n > 1 to find nth root (of value entered above)
Roots by Newton's (iterative) method of successive (improved) approximations
Xn = Xp - f(Xp)/f'(Xp) i.e. Xnext = Xprev - f'(Xprev)/f(Xprev)
.
./|f(Xp)
./ |
./ |
. /___|
. Xn Xp Method to find Xn is derived from:
f(Xp)/(Xp-Xn) = rise/run = slope = f'(Xp)
So to find 4th root of 5,
Use/solve f(x) = x**4 - 5 = 0
'''
print( header )
print()
#import math
def f( x, n, v ):
try:
return x**n - v ## example 4th root of 5 ==> solve: x**4 - 5 = 0 ##
except:
print( "Error in f..." )
return 0
def df( x, n ):
try:
return n*x**(n-1)
except:
print( "Error in df..." )
return 1
def getx1( x, n, v ):
return x - f(x,n,v)/df(x,n)
def more():
return input("More (y/n) ? " ).lower() != 'n'
def takeIn():
done = False
while not done:
floatErr = True
overflowErr = True
try:
v = eval(input( "Enter number to find nth root : " ) )
if v < 0:
raise
floatErr = False
n = int(input( "Enter n for nth root to find : " ) )
if n < 2:
overflowErr = False
raise
case = 0
if v < 1:
if n*v >= 1:
test = v - 1/n - (v - 1/n)**n - v
case = 1
else:
test = v*n - (v*n)**n - v
case = 2
elif v > 1:
test = (1+v/n) - (1+v/n)**n - v
case = 3
overflowErr = False
done = True
except:
if floatErr:
print( "Entry error! Positive number expected ... Try again ... " )
elif overflowErr:
if case == 1: x = v - 1/n
elif case == 2: x = n*v
else: x = 1 + v/n
print( "Overflow error! {}**{} causes overflow ... Try again ... ".format(x,n) )
else:
print( "Entry error! Positive integer > 1 expected ... Try again ... " )
return v, n
while 1:
v, n = takeIn()
if v < 1:
if n*v >= 1:
x1 = v - 1/n # (n*v-1)/n
else:
x1 = n*v
else:
x1 = 1 + v/n
x0 = 0
count = 0
ok = True
while ok:
count += 1
x0 = x1
try:
x1 = getx1(x1,n,v)
#print( count, x0, x1 )
except:
print( "Some error occured ... perhaps overflow ...")
ok = False
if (1 - err) <= x0/x1 <= (1 + err):
break
if count > 10000:
print( "Loops > 10,000 ERROR!" )
break
if n > 3: th = 'th'
elif n == 3: th = 'rd'
else: th = 'nd' # since n == 2 here #
print( '{}{}_root({}) = {} found in {} loops.'.\
format(n, th, v, x1, count) )
print( "Check if f({}) = {} = 0.0 (or nearly so?)".\
format( x1, f(x1,n,v) ) )
if not more():
break
print()
input("Press 'Enter' to continue/exit ... ")