Author Topic: Demo Programs linked from Beginning Calculus  (Read 26593 times)

Offline David

  • Hero Member
  • *****
  • Posts: 644
    • View Profile
Demo Programs linked from Beginning Calculus
« 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


Code: [Select]
# 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')

Code: [Select]
# 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 ...

Code: [Select]
# 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 ... " )
« Last Edit: January 24, 2013, 06:17:41 AM by David »

Offline David

  • Hero Member
  • *****
  • Posts: 644
    • View Profile
Re: Demo Programs linked from Beginning Calculus
« Reply #1 on: January 14, 2013, 06:19:09 AM »
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 )

Code: [Select]
# 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)

Code: [Select]
# 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 ...
'''
« Last Edit: January 18, 2013, 11:13:16 PM by David »

Offline David

  • Hero Member
  • *****
  • Posts: 644
    • View Profile
Re: Demo Programs linked from Beginning Calculus
« Reply #2 on: January 14, 2013, 06:40:02 AM »
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


Code: [Select]
# 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 ...

Code: [Select]
# 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



Offline David

  • Hero Member
  • *****
  • Posts: 644
    • View Profile
Re: Demo Programs linked from Beginning Calculus
« Reply #3 on: January 18, 2013, 11:19:23 PM »
A little demo to show that, as x increases, then e^x exceeds x^n for all finite n ...

Code: [Select]
# 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 ...

Code: [Select]
# 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')
« Last Edit: January 23, 2013, 03:45:17 AM by David »

Offline David

  • Hero Member
  • *****
  • Posts: 644
    • View Profile
Re: Demo Programs linked from Beginning Calculus
« Reply #4 on: January 23, 2013, 03:46:27 AM »
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)


Code: [Select]
# 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 ...

Code: [Select]
# 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 ...

Code: [Select]
# 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) ...

Code: [Select]
# 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 ... ")
« Last Edit: February 02, 2013, 04:48:45 AM by David »