SageWeek07.sagews

Q1.

f(x) = x ^ 3 + 3 * x + sin(x)
f.variables()

(x,)

f(x, y) = x * y + x ^ 2 + y ^ 2
f.variables()

(x, y)

Q2.

f(x) = 3 * x ^ 2 / (x ^ 3 + x - 1)
plot(f, 1, 10)

limit(f,x=1)

x |--> 3

limit(f,x=oo)

x |--> 0

f.limit(x=1)

x |--> 3

f.limit(x=oo)

x |--> 0

bool(f(1) == f.limit(x=1))  # This shows that you can convert a function expression to a boolean.

True

Q3.

plot(f)

p = plot(f, x, 0.8, 10)  # This creates a first plot object
print type(p)

p += plot(f, x, -10, .6)  # This add another plot to the plot (the other side of the discontinuity)
print type(p)
p.show()  # This shows the plot

<class 'sage.plot.graphics.Graphics'> <class 'sage.plot.graphics.Graphics'>

f.denominator()  # This gives the denominator of f

x |--> x^3 + x - 1

rts = f.denominator().roots()  # This finds the roots of the denominator
rts

[(-1/2*(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3)*(I*sqrt(3) + 1) + 1/6*(-I*sqrt(3) + 1)/(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3), 1), (-1/2*(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3)*(-I*sqrt(3) + 1) + 1/6*(I*sqrt(3) + 1)/(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3), 1), ((1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3) - 1/3/(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3), 1)]

realrts = [r[0] for r in rts if r[0].n().is_real()]  # This obtains all the real roots
realrts

[(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3) - 1/3/(1/18*sqrt(31)*sqrt(3) + 1/2)^(1/3)]

limit(f, x=realrts[0], dir="plus")

x |--> +Infinity

limit(f, x=realrts[0], dir="minus")

x |--> -Infinity

limit(f, x=realrts[0])  # This simply shows that no limit exists.

x |--> Infinity

Q4.

f(x) = exp(x)
g(x) = sin(x)
var('a')
L1 = limit(f(x) + g(x), x = a)
L2 = limit(f(x), x = a) + limit(g(x), x = a)
bool(L1 == L2)

a True

f(x) = exp(x)
g(x) = sin(x)
var('a')
L1 = limit(f(x) * g(x), x = a)
L2 = limit(f(x), x = a) * limit(g(x), x = a)
bool(L1 == L2)

a True

f(x) = exp(x)
g(x) = sin(x)
var('a')
L1 = limit(f(x) / g(x), x = a)
L2 = limit(f(x), x = a) / limit(g(x), x = a)
bool(L1 == L2)

a True

Q5. 

data = [[x, abs(sin(x) - x)] for x in srange(0,.05,.00005)]  # Obtains the data as doublets
len(data)  # verify that have required number of points

1000

list_plot(data)  # This is the basic plot
a=list_plot(data, axes_labels=["x", "abs(sin(x)-x)"], color='red')  # This gives some options
show(a)
# This plot seems to indicate that the limit is 1

f(x) = sin(x) / x
f.limit(x = 0)

x |--> 1

Q6.

e.n()

2.71828182845905

data = [[x, (1 + x) ^ (1 / x)] for x in srange(0,.05,.00005)]  # Obtains the data as doublets
len(data)  # verify that have required number of points

1000

p = list_plot(data, legend_label='$(1+x)^{(1/x)}$')  # Including the legend label option
p += list_plot([[x[0], e] for x in data], color='red', legend_label = 'e')
p.show()

Q7.

n = var('n')
f(x) = x ^ n
diff(f,x)

x |--> n*x^(n - 1)

Q8.

f(x) = x ^ 3 + 3 * x - 20

var('h')
expr = ((f(x + h) - f(x)) / h)
bool(expr.expand().limit(h=0) == diff(f,x))  # Students might like to carry this out in steps

h True

Q9.

def tangplot(f, a, x1, x2):
    """
    A function to plot f(x) as well as tangent at x=a

    Arguments:
        f: A function
        a: A point
        x1: Lower bound for domain of plot # This wasn't asked for in this question but would be great if some students realised it would make the plot look nicer
        x2: Upper bound for domain of plot # Same comment as above

    Output: A plot of f(x) as well as the tangent line at x=a
    """
    p = plot(f, x, x1, x2, legend_label = "$f(x)$")
    fdash(x) = diff(f,x)  #  Define a new function for the derivative
    p += plot(fdash(a) * (x - a) + f(a), x, x1, x2, color='red', legend_label = "Tangent at $x=%s$" % a)  # This makes use of some basic algebra
    return p

f(x) = sin(x) + 3 * x + 1 / x
p = tangplot(f, 2, .05, 10)
p

Q10.

f(x) = exp(x)
g(x) = sin(x)
D1 = diff(f(x) + g(x), x)
D2 = diff(f(x), x) + diff(g(x), x)
bool(D1 == D2)

True

f(x) = exp(x)
g(x) = sin(x)
D1 = diff(f(x) * g(x), x)
D2 = diff(f(x), x) * g(x) + f(x) * diff(g(x), x)
bool(D1 == D2)

True

f(x) = exp(x)
g(x) = sin(x)
D1 = diff(f(x) / g(x), x)
D2 = (diff(f(x), x) * g(x) - f(x) * diff(g(x), x)) / (g(x) ^ 2)
bool(D1 == D2)

True

Q11.

f(x) = x ^ 4
integrate(f, x)

x |--> 1/5*x^5

f.integrate(x)

x |--> 1/5*x^5

integrate(f, x, 5, 12)

245707/5

Q12.

a = var('a')
f(x) = sin(a * x)
f.integrate(x)

x |--> -cos(a*x)/a

n = var('n')
assume(n != -1)  # Need this to give required result.
f(x) = x ^ n
f.integrate(x)

x |--> x^(n + 1)/(n + 1)

forget()
assume(n == -1)
f(x) = x ^ n
f.integrate(x)

x |--> log(x)

Q13.

def riemann(f, deltax, a, b):
    """
    A function to return the Riemann integral of a function.

    Arguments:
        f: the function to be integrated
        deltax: the step size
        a: the lower boundary
        b: the upper boundary

    Outputs: The Riemann integral
    """
    xpoints = srange(a, b, deltax)  # Get the collection of points
    r = 0
    for x in xpoints:
        r += deltax * f(x)  # Add area of corresponding rectangle
    return r

f(x) = sin(x) + x ^ 4 - 10
a = 10
b = 11
exact = f.integral(x, a, b)
approxlist = [[deltax, abs(riemann(f, deltax, a,b) - exact)] for deltax in srange(.001, 1, .001)]
list_plot(approxlist, axes_labels=["$\Delta x$", "Difference"])

Q14.

import csv

f = open('W07_D01.csv', 'r')  # Importing the data as shown in the video
data = csv.reader(f)
data = [row for row in data]
f.close()

data = [[eval(i) for i in row] for row in data[1:]]  # Remove first row

class dataintegral():  # There is no need to create a class to do this however using classes is to be encouraged and allows for nice manipulation further on
    """
    A class for each integral.

    Methods: init

    Attributes:
        a: coefficient of x ^ 2
        b: coefficient of x
        c: coefficient of x ^ 0
        A: lower boundary
        B: upper boundary
        integral: the value of the integral
        nonneg: a boolean - True if integral is non neg

    """
    def __init__(self, a, b, c, A, B):
        """
        Initialisation method.

        Arguments:
            a: coefficient of x ^ 2
            b: coefficient of x
            c: coefficient of x ^ 0
            A: lower boundary
            B: upper boundary

        """
        self.a = a
        self.b = b
        self.c = c
        self.A = A
        self.B = B
        self.integral = integrate(a * x ^ 2 + b * x + c, x, self.A, self.B)
        self.nonneg = True
        if self.integral < 0:
            self.nonneg = False

data = [dataintegral(r[0], r[1], r[2], r[3], r[4]) for r in data]

mean([i.integral for i in data]).n()

-43836.9658333333

mean([i.integral for i in data if i.nonneg]).n()

701685.174382716

mean([i.integral for i in data if not i.nonneg]).n()

-919015.130434783