14.1 _ Class ID 26545 _ 1014789.ipynb

Kernel: SageMath 8.4

1 High Dimensional Chains

# a)
x(t) = (t+1)^2
y(t) = 3*t^2

show(x, y)

parametric_plot((x(t), y(t)), (t, 0, 5), 
               axes_labels=['$x$', '$y$'])

T(x,y) = e^(5*x)*(x^2 + y^2)
show(T)

x(t) = (t+1)^2
y(t) = 3*t^2
T_t(t) = e^(5*x)*(x^2 + y^2)
show(T_t)
para_plt = parametric_plot3d((T_t,x, y), (t, 0, 5), boundary_style={"color": "black", "thickness": 2})
para_plt

print 'Lowest value of temperature is %.4f' % T_t(0)

Lowest value of temperature is 148.4132

1.d) Testing chain rule for T(t)

var('t')

# Calculating dT/dt with chain rule
T(x, y) = e^(5*x)*(x^2+y^2)
print 'Partial derivatives of T(t) in terms of x and y'
show(T_1.diff())

dT_dx = T.diff()[0]
dT_dy = T.diff()[1]
print'dT/dt with chain rule'
check_dT_dt = dT_dx((t+1)^2,3*t^2)*dx_dt + dT_dy((t+1)^2,3*t^2)*dy_dt
show(check_dT_dt)


# Calculating dT/dt without chain rule
x(t) = (t+1)^2
y(t) = 3*t^2

# Recreating T(t) in terms of t only.
T(t) = e^(5*x)*(x^2+y^2)
print 'T(t) in terms of t'
show(T(t))
test_dT_dt = T.diff()
print'dT/dt without chain rule'
show(test_dT_dt)

# Equating those two equations
print 'Checking above two equations'
solve(check_dT_dt - test_dT_dt == 0, t)

Partial derivatives of T(t) in terms of x and y

dT/dt with chain rule

T(t) in terms of t

dT/dt without chain rule

Checking above two equations

[t == t]

2. Taking up a notch

# Plotting regtion E{(x, y, z)|0 <=z <=2 - x^2/2, 0 <=y <= x}
var('x, y, z')

region1 = implicit_plot3d(y==x, (x, 0, 2), (y, 0, 2), (z, 0, 2), 
                         plot_points=100, 
                         color = 'red',
                          legend_label = '$y<=x$',
                          axes_labels=['$x$','$y$','$z$'],
                         region=lambda x,y,z:  x>=0 and x<=2 and z>= 0 and z<=2-x^2/2)

region2 = implicit_plot3d(z==2-x^2/2, (x, 0, 2), (y, 0, 2), (z, 0, 2), 
                         plot_points=100, 
                          legend_label = '$z<=2-x^2/2$',
                         region=lambda x,y,z:  y <= x and y >=0 and x>=0 and x<=2)
region1 + region2

# triple integral of 1
var('x, y, z')
f(x, y, z) = 1 
# Triple integration over f(x, y, z) = 1 for calculating volme
v = integral(integral(integral(f(x, y ,z), z, 0, 2-x^2/2), y, 0, x), x, 0, 2)
print 'Volume or Triple Integral of f(x, y, z) = 1 over E is %.2f' % v

# Calculating mass with triple integrals over density function
var('x, y, z')
var('k')
f(x, y, z) = k * z
m = integral(integral(integral(f(x, y, z), z, 0, 2-x^2/2), y, 0, x), x, 0, 2)
print 'Mass or Triple Integral of f(x, y, z) = kz over E is %s' % m

Volume or Triple Integral of f(x, y, z) = 1 over E is 2.00
Mass or Triple Integral of f(x, y, z) = kz over E is 4/3*k

# Calculating Center of mass
var('x, y, z')
var('k')
d(x, y, z) = k*z

m_yz = integral(integral(integral(x*d, z, 0, 2-x^2/2), y, 0, x), x, 0, 2)
print 'M_yz: %s' % m_yz

m_xz = integral(integral(integral(y*d, z, 0, 2-x^2/2), y, 0, x), x, 0, 2)
print 'M_xz: %s' % m_xz

m_xy = integral(integral(integral(z*d, z, 0, 2-x^2/2), y, 0, x), x, 0, 2)
print 'M_xy: %s' % m_xy

x_center = m_yz/m
y_center = m_xz/m
z_center = m_xy/m

print 'The center of mass is: (%s, %s, %s)' % (x_center, y_center, z_center)

M_yz: 128/105*k
M_xz: 64/105*k
M_xy: 4/3*k
The center of mass is: (32/35, 16/35, 1)

3 Euler Solver

# Algorithm for solving estimated value of y

def savings_estimator(r=0.01, q=250, 
                      n = 12, y0=100, h=1, 
                      t0=0, verbose=False, 
                      return_arr=False):
    var('y0, t0, h')

    #y1 = y0 + f(t0, y0)*h
    '''
        Inputs:
            r -> Interest rate (0.01 in question)
            q -> savings rate (250$ per month in question)
            y0 -> Initial investment (Default y = 100)
            t0 -> Initial time (Default t=0)
            h -> step size
            n -> number of months (Default n=12 or one year)
            verbose -> Denotes whether or not to print savings after each
                        step size
            return_arr -> Status to return the arrays storing savings and
                        time in each increment
        Outputs:
            - Prints the savings after each increment
            - Returns the savings after a year
            
    '''

    y_s = [y0] # array to store savings after each incresment
    t_s = [t0] # array to store time steps
    
    f(t, y) = r*y + q #dy/dt equation

    i = 1 # index position of array
    t_n = t0 # temparay variable to store time steps

    while int(t_n) < n:
        #t_n = t_n-1 + h
        t_n = t_s[i-1] + h
        
        #y_n = y_n-1 + f(t_n-1, y_n-1)*h
        y_n = y_s[i-1] + f(t_s[i-1], y_s[i-1])*h
        y_s.append(y_n)
        t_s.append(t_n)

        if verbose:
            # Prints savings after each step size
            print 'Amount at t=%.2f is %.2f' % (t_s[i], y_s[i])
        
        i += 1
        
    if return_arr:
        # Returns all the arrays if required
        return t_s, y_s
    
    return y_s[-1]

# Using following values of step sizes.
h_arr = [1, 2, 3, 0.1, 0.5, 0.01, 0.001, 0.0005]

h_arr.sort(reverse=True)

for h in h_arr:
    print 'For h=%f, savings after a year = %.2f' % (h, savings_estimator(h=h))

For h=3.000000, savings after a year = 3250.27
For h=2.000000, savings after a year = 3266.68
For h=1.000000, savings after a year = 3283.31
For h=0.500000, savings after a year = 3291.71
For h=0.100000, savings after a year = 3326.77
For h=0.010000, savings after a year = 3302.83
For h=0.001000, savings after a year = 3300.44
For h=0.000500, savings after a year = 3300.16

# plot slope field

r = 0.01 # interest rate
q = 250 # savings rate

dy_dt = r*y + q
plt_slope = plot_slope_field(dy_dt, (t,0,12), (y,0,3500),
                            axes_labels=['$t$ (in months)',
                                         '$y$-Savings in dollars'])

# getting output from savings estimator
step_size = 1
print 'Step_size: %.2f' % step_size
(t_s, y_s) = savings_estimator(h=step_size, verbose=True, return_arr=True)

# storing the time and savings in each step to points array
points  = []
for ind in range(len(t_s)):
    points.append([t_s[ind], y_s[ind]])

plt_point = scatter_plot(points)

show(plt_slope + plt_point)

Step_size: 1.00
Amount at t=1.00 is 351.00
Amount at t=2.00 is 604.51
Amount at t=3.00 is 860.56
Amount at t=4.00 is 1119.16
Amount at t=5.00 is 1380.35
Amount at t=6.00 is 1644.16
Amount at t=7.00 is 1910.60
Amount at t=8.00 is 2179.70
Amount at t=9.00 is 2451.50
Amount at t=10.00 is 2726.02
Amount at t=11.00 is 3003.28
Amount at t=12.00 is 3283.31

# Due to memory and implementation issues this code does not run in CoCalc
# Please run it on the local SageMath server of your laptop
t = var('t')
y = function('y')(t)
r = 0.01 # interest rate
q(t) = 250 # monthly savings rate

soln = desolve(diff(y, t) -r*y - q(t), y)

print 'Solution of differential equation is'
show(soln)
# Solution appears as (C−25000*e(−t/100))e^(t/100)

print 'In simpler terms: y(t) is'
var('C')
y_actual(t) = C*e^(t/100) - 25000
show(y_actual(t))

# Finding value of C using intial values
# y(0) = 100, when t = 0
results = solve(y_actual(0) == 100, C)
C = results[0].rhs()
print 'The value of C is %.2f' % C

# Redefining the function with value of C
y_actual(t) = C*e^(t/100) - 25000
print 'The actual function y(t) is '
show(y_actual(t))

print 'Actual savings after a year is: %.4f' % y_actual(12)

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-65-37366c582a9b> in <module>()
      6 __tmp__=var("t"); q = symbolic_expression(Integer(250) ).function(t)# monthly savings rate
      7 
----> 8 soln = desolve(diff(y, t) -r*y - q(t), y)
      9 
     10 print 'Solution of differential equation is'
/ext/sage/sage-8.4_1804/local/lib/python2.7/site-packages/sage/calculus/desolvers.pyc in desolve(de, dvar, ics, ivar, show_method, contrib_ode, algorithm)
    595                 raise NotImplementedError("Maxima was unable to solve this ODE.")
    596         else:
--> 597             raise NotImplementedError("Maxima was unable to solve this ODE. Consider to set option contrib_ode to True.")
    598 
    599     if show_method:
NotImplementedError: Maxima was unable to solve this ODE. Consider to set option contrib_ode to True.

4 Something fishy

# Plotting slope field for dP/dt = rP(1-P/N)
t, P = var('t,P')
N = 100 # max. capacity of 10,000 fish
r = 0.14 # reproductive rate of 14% per week
plt_slope = plot_slope_field(r*P*(1-(P/N)),
                             (t,0,52),(P,-2,N+2),
                             axes_labels=['$t$', '$P$']
                            )
plt_line = plot(N,(x,0,52))

plt_slope + plt_line

t, P, H = var('t,P, H')

N = 100 # max. capacity of 10,000 fish
r = 0.14 # reproductive rate of 14% per week
plt_slope = plot_slope_field((r*P*(1-(P/N))-H),
                             (H,0,5), (P,0,N+2),
                            axes_labels=['$H$', '$P$'])

plt_slope.show()

plt_slope = plot_slope_field((r*P*(1-(P/N))-H),
                             (P,0,N), (H,0,5), 
                            axes_labels=['$P$', '$H$'])

plt_slope.show()


H = 4
plt_slope = plot_slope_field((r*P*(1-(P/N))-H),
                             (t,0,100),(P,0,N),
                            axes_labels=['$t$', '$P$'])

plt_slope.show()

var('N, H, r')
N=100
r=0.14
P_plus(H) = (N + sqrt(N^2 - 4*H*N/r))/2
P_minus(H) = (N - sqrt(N^2 - 4*H*N/r))/2

print P_plus(3)
print P_minus(3)

68.8982236504614
31.1017763495386

Taylor Approximation

# Linear approximation
f(x) = sin(x)

x_0 = 0
lin_approx(x) = f(x_0) + f.diff()(x)*(x - x_0)

plot_f = plot(f(x), (x, -1, 1), 
              axes_labels=['$x$','$y$'], 
              color='blue', legend_label='sin(x)')
plot_lin = plot(lin_approx(x), (x, -1, 1), color='red',
               legend_label = 'Linear approximation')

# plotting sin(x) and its linear approximations at x = 0,
show(plot_f + plot_lin)



# Calculating error
print 'sin(pi/6) = %.4f' % (sin(pi/6))
print 'Linear approximation at x = pi/6 = %.4f' % (lin_approx(pi/6))

print 'Absolute error = %.4f' % abs(sin(pi/6)-lin_approx(pi/6))

# quadratic approximation

var('a, b, c, x')

p(x) = a*x^2 + b*x + c
p.diff()

p_0 = p(0)
print 'P(0)='
show(p_0)

p_1st = p.diff()
print 'P_lst(0) = '
show(p_1st(0))

p_2nd  = p_1st.diff()
print 'P_2nd(0) = '
show(p_2nd(0))
# thus a = p_2nd(0)/2


# Main equation
f(x) = sin(x)

f_1st = f.diff()

f_2nd = f_1st.diff()


f_quad(x) = (f_2nd(0)/2)*x^2 + f_1st(0)*x + f(0)

print 'Quadratic approximation'
show(f_quad(x))

# Ploting sin(x) and its quadriatic approximation

plot_sin = plot(sin(x), (x, -1, 1), 
                axes_labels=['$x$', '$y$'], 
                legend_label='$f(x) = sin(x)$', 
                color='blue'
               )

plot_quad = plot(f_quad(x), (x, -1, 1),
                 axes_labels=['$x$', '$y$'], 
                 legend_label='$Quad approx = x$', 
                 color='red'
                )
show(plot_sin + plot_quad)


# Error calculation
print 'sin(pi/6) = %.4f' % (f(pi/6))
print 'Quadratic approximation of sin(pi/6) = %.4f' % (f_quad(pi/6))
print 'Error = %.4f' % abs(f_quad(pi/6) - f(pi/6))

# Cubic approximation
var('a, b, c, d')
p(x) = a*x^3 + b*x^2+ c*x + d
print 'p(0) = %s' % p(0)
p_1st = p.diff()
print 'p`(0) = %s' % p_1st(0)
p_2nd = p_1st.diff()
print 'p``(0) = %s' % p_2nd(0)
# b = p_2nd(0)/2
p_3rd = derivative(p_2nd, x)
print 'p```(0) = %s' % p_3rd(0)
# a = p_3rd(0)/6


f(x) = sin(x)
f_1st = derivative(f(x), x)
f_2nd = derivative(f_1st, x)
f_3rd = derivative(f_2nd, x)

d = f(0)
c = f_1st(0)
b = f_2nd(0)/2
a = f_3rd(0)/6

f_cubic(x) = a*x^3 + b*x^2 + c*x + d
show(f_cubic)

# Ploting sin(x) and its cubic approximation

plot_sin = plot(f(x), (x, -2, 2), 
                axes_labels=['$x$', '$y$'], 
                legend_label='$f(x) = sin(x)$', 
                color='blue'
               )

plot_cubic = plot(f_cubic(x), (x, -2, 2),
                 axes_labels=['$x$', '$y$'], 
                 legend_label='$Cubic approximation$', 
                 color='red'
                )
show(plot_sin + plot_cubic)


# Error calculation
print 'sin(pi/6) = %.4f' % (f(pi/6))
print 'Quadratic approximation of sin(pi/6) = %.4f' % (f_cubic(pi/6))
print 'Error = %.4f' % abs(f_cubic(pi/6) - f(pi/6))

# Taylor polynomials
# g)

def plot_taylor(func, x_range):
    f(x) = func
    x_range = x_range
    plot_all = plot(f(x), x_range, 
                    axes_labels=['$x$', '$y$'], 
                    legend_label=func, 
                    color='blue'
                   )

    color_taylor = ['orange',
                    'red',
                    'pink',
                    'green',
                    'yellow',
                    'black'
                   ]

    for i in range(1, 7):
        color = color_taylor[i-1]
        legend = '%s approximation' % i
        print '%d degree' % i
        show(f.taylor(x, 0, i))
        plot_i = plot(f.taylor(x, 0, i), x_range,
                      axes_labels=['$x$', '$y$'], 
                      legend_label=legend, 
                      color=color
                     )
        plot_all += plot_i


    show(plot_all)

print 'Taylor approximation of f(x)=sin(x)'
f(x) = sin(x)
plot_taylor(f(x), (x, -2*pi, 2*pi))
print '--'*20

print 'Taylor approximation of g(x)=e^x'
g(x) = e^x
plot_taylor(g(x), (x, -2*pi, 2*pi))
print '--'*20

print 'Taylor approximation of h(x) = 1/(1-x)'
h(x) = 1/(1-x)
plot_taylor(h(x), (x, -5, 5))
print '--'*20

Taylor approximation of f(x)=sin(x)
1 degree

2 degree

3 degree

4 degree

5 degree

6 degree

----------------------------------------
Taylor approximation of g(x)=e^x
1 degree

2 degree

3 degree

4 degree

5 degree

6 degree

----------------------------------------
Taylor approximation of h(x) = 1/(1-x)
1 degree

2 degree

3 degree

4 degree

5 degree

6 degree

----------------------------------------

# (i)

print 'f(x)'
f(x) = (1/(sqrt(2*pi)))*e^(-x^2/2)
show(f(x))

print '1st Taylor Polynomials'
p_1 = f.taylor(x, 0, 1)
show(p_1)
print 'Integrate p_1 from [-1, 1] = %.5f ' % (integrate(p_1, x, -1, 1))
show(integrate(p_1, x, -1, 1))

print '2nd Taylor Polynomials'
p_2 = f.taylor(x, 0, 2)
show(p_2)
print 'Integrate p_2 from [-1, 1] = %.5f ' % (integrate(p_2, x, -1, 1))
show(integrate(p_2, x, -1, 1))

print '3rd Taylor Polynomials'
p_3 = f.taylor(x, 0, 3)
show(p_3)
print 'Integrate p_3 from [-1, 1] = %.5f ' % (integrate(p_3, x, -1, 1))
show(integrate(p_3, x, -1, 1))

f(x)

1st Taylor Polynomials

Integrate p_1 from [-1, 1] = 0.79788 

2nd Taylor Polynomials

Integrate p_2 from [-1, 1] = 0.66490 

3rd Taylor Polynomials

Integrate p_3 from [-1, 1] = 0.66490 

f(x) = (1/(sqrt(2*pi)))*e^(-x^2/2)
show(f(x))
for i in range(1, 20):
#     print '%d degree Taylor Polynomials' % i
    f_approx = f.taylor(x, 0, i)
    print 'Integration %d degree Taylor Polynomial from [-1, 1] = %.4f ' % (i, integrate(f_approx, x, -1, 1))


print '--'*20
print 'Integration of f(x) [-1, 1] = %.4f ' % (integrate(f(x), x, -1, 1))

Integrate 1 degree Taylor Polynomial from [-1, 1] = 0.7979 
Integrate 2 degree Taylor Polynomial from [-1, 1] = 0.6649 
Integrate 3 degree Taylor Polynomial from [-1, 1] = 0.6649 
Integrate 4 degree Taylor Polynomial from [-1, 1] = 0.6849 
Integrate 5 degree Taylor Polynomial from [-1, 1] = 0.6849 
Integrate 6 degree Taylor Polynomial from [-1, 1] = 0.6825 
Integrate 7 degree Taylor Polynomial from [-1, 1] = 0.6825 
Integrate 8 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 9 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 10 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 11 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 12 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 13 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 14 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 15 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 16 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 17 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 18 degree Taylor Polynomial from [-1, 1] = 0.6827 
Integrate 19 degree Taylor Polynomial from [-1, 1] = 0.6827 
----------------------------------------
Integrate of f(x) [-1, 1] = 0.6827 

f(x) = (1/(sqrt(2*pi)))*e^(-x^2/2)

f_approx = f.taylor(x, 0, 10)
print 'Integration of f(x) [-1, 1] = %.4f ' % (integrate(f(x), x, -oo, oo))
# print 'Integrate %d degree Taylor Polynomial from [-1, 1] = %.4f ' % (10, integrate(f_approx, x, -oo, oo))

show(f_approx)

Integration of f(x) [-1, 1] = 1.0000 

f_first(x) = x
f_third(x)= -x^3/6+x
f_fifth(x)= x^5/120 -x^3/6+x

# Error calculation for pi/6
print 'sin(pi/6) = %.4f' % (f(pi/6))
print 'First and second order approximation of sin(pi/6) = %.4f' % (f_first(pi/6))
print 'Error = %.4f' % abs(f_first(pi/6) - f(pi/6))

print 'sin(pi/6) = %.4f' % (f(pi/6))
print 'Third and fourth order approximation of sin(pi/6) = %.4f' % (f_third(pi/6))
print 'Error = %.4f' % abs(f_third(pi/6) - f(pi/6))

print 'sin(pi/6) = %.4f' % (f(pi/6))
print 'Fifth and sixth order approximation of sin(pi/6) = %.4f' % (f_fifth(pi/6))
print 'Error = %.4f' % abs(f_fifth(pi/6) - f(pi/6))

# Error calculation for 8pi

print 'sin(8pi) = %.4f' % (f(8*pi))
print 'First and second order approximation of sin(8pi) = %.4f' % (f_first(8*pi))
print 'Error = %.4f' % abs(f_first(8*pi) - f(8*pi))

print 'sin(8pi) = %.4f' % (f(8*pi))
print 'Third and fourth order approximation of sin(8pi) = %.4f' % (f_third(8*pi))
print 'Error = %.4f' % abs(f_third(8*pi) - f(8*pi))

print 'sin(8pi) = %.4f' % (f(8*pi))
print 'Fifth and sixth order approximation of sin(8*pi) = %.4f' % (f_fifth(8*pi))
print 'Error = %.4f' % abs(f_fifth(8*pi) - f(8*pi))

sin(pi/6) = 0.5000
First and second order approximation of sin(pi/6) = 0.5236
Error = 0.0236
sin(pi/6) = 0.5000
Third and fourth order approximation of sin(pi/6) = 0.4997
Error = 0.0003
sin(pi/6) = 0.5000
Fifth and sixth order approximation of sin(pi/6) = 0.5000
Error = 0.0000
sin(8pi) = 0.0000
First and second order approximation of sin(8pi) = 25.1327
Error = 25.1327
sin(8pi) = 0.0000
Third and fourth order approximation of sin(8pi) = -2620.7362
Error = 2620.7362
sin(8pi) = 0.0000
Fifth and sixth order approximation of sin(8*pi) = 80943.0391
Error = 80943.0391

f_first(x) = x+1
f_second(x) = (1/2)*x^2+x+1
f_third(x)= (1/6)*x^3+(1/2)*x^2+x+1
f_fourth(x) = (1/24)*x^4+(1/6)*x^3+(1/2)*x^2+x+1
f_fifth(x)= x^5/120+(1/24)*x^4+(1/6)*x^3+(1/2)*x^2+x+1
f_sixth(x) = (1/720)*x^6+x^5/120+(1/24)*x^4+(1/6)*x^3+(1/2)*x^2+x+1

# Error calculation for 0.01
print 'e^(0.01) = %.4f' % (e^(0.01))
print 'First order approximation of e^(0.01) = %.4f' % (f_first(0.01))
print 'Error = %.4f' % abs(f_first(0.01) - e^(0.01))

print 'e^(0.01) = %.4f' % (e^(0.01))
print 'Second order approximation of e^(0.01) = %.4f' % (f_second(0.01))
print 'Error = %.4f' % abs(f_second(0.01) - e^(0.01))

print 'e^(0.01) = %.4f' % (e^(0.01))
print 'Third order approximation of e^(0.01) = %.4f' % (f_third(0.01))
print 'Error = %.4f' % abs(f_third(0.01) - e^(0.01))

print 'e^(0.01) = %.4f' % (e^(0.01))
print 'Fourth order approximation of e^(0.01) = %.4f' % (f_fourth(0.01))
print 'Error = %.4f' % abs(f_fourth(0.01) - e^(0.01))

print 'e^(0.01) = %.4f' % (e^(0.01))
print 'Fifth order approximation of e^(0.01) = %.4f' % (f_fifth(0.01))
print 'Error = %.4f' % abs(f_fifth(0.01) - e^(0.01))

print 'e^(0.01) = %.4f' % (e^(0.01))
print 'Sixth order approximation of e^(0.01) = %.4f' % (f_sixth(0.01))
print 'Error = %.4f' % abs(f_sixth(0.01) - e^(0.01))

# Error calculation for 100

print 'e^(100) = %.4f' % (e^(100))
print 'First order approximation of e^(1000) = %.4f' % (f_first(100))
print 'Error = %.4f' % abs(f_first(100) - e^(100))

print 'e^(100) = %.4f' % (e^(100))
print 'Second order approximation of e^(100) = %.4f' % (f_second(100))
print 'Error = %.4f' % abs(f_second(100) - e^(100))

print 'e^(100) = %.4f' % (e^(100))
print 'Third order approximation of e^(100) = %.4f' % (f_third(100))
print 'Error = %.4f' % abs(f_third(100) - e^(100))

print 'e^(100) = %.4f' % (e^(100))
print 'Fourth order approximation of e^(100) = %.4f' % (f_fourth(100))
print 'Error = %.4f' % abs(f_fourth(100) - e^(100))

print 'e^(100) = %.4f' % (e^(100))
print 'Fifth order approximation of e^(1000) = %.4f' % (f_fifth(100))
print 'Error = %.4f' % abs(f_fifth(100) - e^(100))

print 'e^(100) = %.4f' % (e^(100))
print 'Sixth order approximation of e^(100) = %.4f' % (f_sixth(100))
print 'Error = %.4f' % abs(f_sixth(100) - e^(100))

e^(0.01) = 1.0101
First order approximation of e^(0.01) = 1.0100
Error = 0.0001
e^(0.01) = 1.0101
Second order approximation of e^(0.01) = 1.0101
Error = 0.0000
e^(0.01) = 1.0101
Third order approximation of e^(0.01) = 1.0101
Error = 0.0000
e^(0.01) = 1.0101
Fourth order approximation of e^(0.01) = 1.0101
Error = 0.0000
e^(0.01) = 1.0101
Fifth order approximation of e^(0.01) = 1.0101
Error = 0.0000
e^(0.01) = 1.0101
Sixth order approximation of e^(0.01) = 1.0101
Error = 0.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
First order approximation of e^(1000) = 101.0000
Error = 26881171418161356094253400435962903554686976.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
Second order approximation of e^(100) = 5101.0000
Error = 26881171418161356094253400435962903554686976.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
Third order approximation of e^(100) = 171767.6667
Error = 26881171418161356094253400435962903554686976.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
Fourth order approximation of e^(100) = 4338434.3333
Error = 26881171418161356094253400435962903554686976.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
Fifth order approximation of e^(1000) = 87671767.6667
Error = 26881171418161356094253400435962903554686976.0000
e^(100) = 26881171418161356094253400435962903554686976.0000
Sixth order approximation of e^(100) = 1476560656.5556
Error = 26881171418161356094253400435962903554686976.0000

f_first(x) = x+1
f_second(x) = x^2+x+1
f_third(x)= x^3+x^2+x+1
f_fourth(x) = x^4+x^3+x^2+x+1
f_fifth(x)= x^5+x^4+x^3+x^2+x+1
f_sixth(x) = x^6+x^5+x^4+x^3+x^2+x+1

# Error calculation for 1.01
print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'First order approximation of 1/(1-1.01) = %.4f' % (f_first(1.01))
print 'Error = %.4f' % abs(f_first(1.01) - 1/(1-1.01))

print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'Second order approximation of 1/(1-1.01) = %.4f' % (f_second(1.01))
print 'Error = %.4f' % abs(f_second(1.01) - 1/(1-1.01))

print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'Third order approximation of 1/(1-1.01) = %.4f' % (f_third(1.01))
print 'Error = %.4f' % abs(f_third(1.01) - 1/(1-1.01))

print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'Fourth order approximation of 1/(1-1.01) = %.4f' % (f_fourth(1.01))
print 'Error = %.4f' % abs(f_fourth(1.01) - 1/(1-1.01))

print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'Fifth order approximation of 1/(1-1.01) = %.4f' % (f_fifth(1.01))
print 'Error = %.4f' % abs(f_fifth(1.01) - 1/(1-1.01))

print '1/(1-1.01)) = %.4f' % (1/(1-1.01))
print 'Sixth order approximation of 1/(1-1.01) = %.4f' % (f_sixth(1.01))
print 'Error = %.4f' % abs(f_sixth(1.01) - 1/(1-1.01))


# Error calculation for 0.99

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'First order approximation of 1/(1-0.99) = %.4f' % (f_first(0.99))
print 'Error = %.4f' % abs(f_first(0.99) - 1/(1-0.99))

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'Second order approximation of 1/(1-0.99) = %.4f' % (f_second(0.99))
print 'Error = %.4f' % abs(f_second(0.99) - 1/(1-0.99))

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'Third order approximation of 1/(1-0.99) = %.4f' % (f_third(0.99))
print 'Error = %.4f' % abs(f_third(0.99) - 1/(1-0.99))

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'Fourth order approximation of 1/(1-0.99) = %.4f' % (f_fourth(0.99))
print 'Error = %.4f' % abs(f_fourth(0.99) - 1/(1-0.99))

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'Fifth order approximation of 1/(1-0.99) = %.4f' % (f_fifth(0.99))
print 'Error = %.4f' % abs(f_fifth(0.99) - 1/(1-0.99))

print '1/(1-0.99)) = %.4f' % (1/(1-0.99))
print 'Sixth order approximation of 1/(1-0.99) = %.4f' % (f_sixth(0.99))
print 'Error = %.4f' % abs(f_sixth(0.99) - 1/(1-0.99))

# Error calculation for 100

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'First order approximation of 1/(1-100) = %.4f' % (f_first(100))
print 'Error = %.4f' % abs(f_first(100) - 1/(1-100))

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'Second order approximation of 1/(1-100) = %.4f' % (f_second(100))
print 'Error = %.4f' % abs(f_second(100) - 1/(1-100))

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'Third order approximation of 1/(1-100) = %.4f' % (f_third(100))
print 'Error = %.4f' % abs(f_third(100) - 1/(1-100))

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'Fourth order approximation of 1/(1-100) = %.4f' % (f_fourth(100))
print 'Error = %.4f' % abs(f_fourth(100) - 1/(1-100))

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'Fifth order approximation of 1/(1-100) = %.4f' % (f_fifth(100))
print 'Error = %.4f' % abs(f_fifth(100) - 1/(1-100))

print '1/(1-100)) = %.4f' % (1/(1-100))
print 'Sixth order approximation of 1/(1-100) = %.4f' % (f_sixth(100))
print 'Error = %.4f' % abs(f_sixth(100) - 1/(1-100))

1/(1-1.01)) = -100.0000
First order approximation of 1/(1-1.01) = 2.0100
Error = 102.0100
1/(1-1.01)) = -100.0000
Second order approximation of 1/(1-1.01) = 3.0301
Error = 103.0301
1/(1-1.01)) = -100.0000
Third order approximation of 1/(1-1.01) = 4.0604
Error = 104.0604
1/(1-1.01)) = -100.0000
Fourth order approximation of 1/(1-1.01) = 5.1010
Error = 105.1010
1/(1-1.01)) = -100.0000
Fifth order approximation of 1/(1-1.01) = 6.1520
Error = 106.1520
1/(1-1.01)) = -100.0000
Sixth order approximation of 1/(1-1.01) = 7.2135
Error = 107.2135
1/(1-0.99)) = 100.0000
First order approximation of 1/(1-0.99) = 1.9900
Error = 98.0100
1/(1-0.99)) = 100.0000
Second order approximation of 1/(1-0.99) = 2.9701
Error = 97.0299
1/(1-0.99)) = 100.0000
Third order approximation of 1/(1-0.99) = 3.9404
Error = 96.0596
1/(1-0.99)) = 100.0000
Fourth order approximation of 1/(1-0.99) = 4.9010
Error = 95.0990
1/(1-0.99)) = 100.0000
Fifth order approximation of 1/(1-0.99) = 5.8520
Error = 94.1480
1/(1-0.99)) = 100.0000
Sixth order approximation of 1/(1-0.99) = 6.7935
Error = 93.2065
1/(1-100)) = -0.0101
First order approximation of 1/(1-100) = 101.0000
Error = 101.0101
1/(1-100)) = -0.0101
Second order approximation of 1/(1-100) = 10101.0000
Error = 10101.0101
1/(1-100)) = -0.0101
Third order approximation of 1/(1-100) = 1010101.0000
Error = 1010101.0101
1/(1-100)) = -0.0101
Fourth order approximation of 1/(1-100) = 101010101.0000
Error = 101010101.0101
1/(1-100)) = -0.0101
Fifth order approximation of 1/(1-100) = 10101010101.0000
Error = 10101010101.0101
1/(1-100)) = -0.0101
Sixth order approximation of 1/(1-100) = 1010101010101.0000
Error = 1010101010101.0101

t(y)= 100*ln((y/100)+250)
solve(t==0,y)

[y == -24900]

# Question (i)

def ceil(quot):
    # Rounds off a decimal number to nearest greater number to the function
    return int(quot) + 1

def floor(quot):
    # Rounds off a decimal number to nearest smaller number to the function
    return int(quot)

def sin_calculator(x, f_approx):
    # Divide the number by pi
    quot = float(x/pi)
    
    if quot % 2 == 0:
        # If the quotient is an even number
        # subtract x from the product of pi and floor value of the quotient
        conv_x = x - floor(quot)*pi
    else:
        # If the quotient is an even number
        # subtract x from the product of pi and ceiling value of the quotient
        conv_x = x - ceil(quot)*pi
    
    # Insert the new converted number the linear approximation function, 
    # stored in the calculator
    return f_approx(conv_x)

def test_calc(x_input):
    # Finding approximation function of sin(x)
    var('x')
    f(x) = sin(x)
    f_approx = f.taylor(x, 0, 1)

    print 'With linear approximation, sin(%.2f) = %.6f' % (x_input, sin_calculator(x_input, f_approx))

    f_approx = f.taylor(x_input, 0, 6)
    print 'With sixth degree taylor approximation, sin(%.2f) = %.6f' % (x_input, sin_calculator(x_input, f_approx))

    print 'Actual value of sin(%.2f) = %.6f' % (x_input, sin(x_input))
    
    print '---'*20
    
test_calc(100)
test_calc(200)

With linear approximation, sin(100.00) = -0.530965
With sixth degree taylor approximation, sin(100.00) = -0.506366
Actual value of sin(100.00) = -0.506366
------------------------------------------------------------
With linear approximation, sin(200.00) = -1.061930
With sixth degree taylor approximation, sin(200.00) = -0.873297
Actual value of sin(200.00) = -0.873297
------------------------------------------------------------

print ceil(31.3)
print ceil(31.7)
print floor(32.7)

32
32
32

5%2

1

int(31.3)

31