#!/usr/bin/env python # coding: utf-8 # # Calculus # In[ ]: from sympy import * from IPython.display import display # This, with init_printing(), will let us "print" to MathJax init_printing() # For each exercise, fill in the function according to its docstring. # In[ ]: x, y, z = symbols('x y z') # ## Derivatives # Compute the following # # $$ \frac{d}{dx}\sin(x)e^x$$ # $$ \frac{\partial}{\partial x}\sin(xy)e^x $$ # $$ \frac{\partial^2}{\partial x\partial y}\sin(xy)e^x $$ # # In[ ]: # In[ ]: # In[ ]: # Recall l'Hopital's rule, which states that if $$\lim_{x\to x_0}\frac{f(x)}{g(x)}$$ is $\frac{0}{0}$, $\frac{\infty}{\infty}$, or $-\frac{\infty}{\infty}$, then it is equal to $$\lim_{x\to x_0} \frac{f'(x)}{g'(x)}$$ (we will not consider other indeterminate forms here). # # Write a function that computes $\lim_{x\to x_0}\frac{f(x)}{g(x)}$. Use the `fraction` function to get the numerator and denominator of an expression, for example # In[ ]: fraction(x/y) # You may assume that the only indeterminate forms are the ones mentioned above, and that l'Hopital's rule will terminate after a finite number of steps. Do not use `limit` (use `subs`). Remember that after taking the derivatives, you will need to put the expression into the form $\frac{f(x)}{g(x)}$ before applying l'Hopital's rule again (what function did we learn that does this?). # In[ ]: def lhopital(expr, x, x0): """ Computes limit(expr, x, x0) using l'Hopital's rule. >>> lhopital(sin(x)/x, x, 0) 1 >>> lhopital(exp(x)/x**2, x, oo) oo >>> lhopital((x**2 - 4*x + 4)/(2 - x), x, 2) 0 >>> lhopital(cos(x), x, 0) 1 >>> lhopital((x + sin(x))/x, x, 0) 2 """ # In[ ]: lhopital(sin(x)/x, x, 0) # In[ ]: lhopital(exp(x)/x**2, x, oo) # In[ ]: lhopital((x**2 - 4*x + 4)/(2 - x), x, 2) # In[ ]: lhopital(cos(x), x, 0) # In[ ]: lhopital((x + sin(x))/x, x, 0) # ## Integrals # Recall that the mean value of a function on an interval $[a, b]$ can be computed by the formula # # $$\frac{1}{b - a}\int_{a}^{b} f{\left (x \right )} dx. % Why doesn't \, work? $$ # # Write a function that computes the mean value of an expression on a given interval. # In[ ]: def average_value(expr, x, a, b): """ Computes the average value of expr with respect to x on [a, b]. >>> average_value(sin(x), x, 0, pi) 2/pi >>> average_value(x, x, 2, 4) 3 >>> average_value(x*y, x, 2, 4) 3*y """ # In[ ]: average_value(sin(x), x, 0, pi) # In[ ]: average_value(x, x, 2, 4) # In[ ]: average_value(x*y, x, 2, 4) # Write a function that takes a list of expressions and produces an "integral table", like # # 1. $$\int \sin(x)dx = -\cos(x) + C$$ # 2. $$\int \cos(x)dx = \sin(x) + C$$ # 3. $$\int e^xdx = e^x + C$$ # 4. $$\int \log(x)dx = x(\log(x) - 1) + C$$ # In[ ]: def int_table(exprs, x): """ Produces a nice integral table of the integrals of exprs >>> int_table([sin(x), cos(x), exp(x), log(x)], x) ⌠ ⎮ sin(x) dx = C - cos(x) ⌡ ⌠ ⎮ cos(x) dx = C + sin(x) ⌡ ⌠ ⎮ x x ⎮ ℯ dx = C + ℯ ⌡ ⌠ ⎮ log(x) dx = C + x⋅log(x) - x ⌡ """ # In[ ]: int_table([sin(x), cos(x), exp(x), log(x)], x) # Now use your function to compute the integrals in this Mathematica ad. Remember that the inverse trig functions are spelled like `asin` in SymPy. # # The ad below probably has a typo, because one of the integrals is trivial to compute. Include what you think the integral should be, and see if SymPy can compute that as well. # In[ ]: from IPython.core.display import Image Image(filename='../imgs/Mathematica-ring-a.png') # In[ ]: # ## Limits # Recall that the definition of the derivative of $f(x)$ at $x=x_0$ is $$f'(x_0) = \lim_{x\to x_0}\frac{f(x) - f(x_0)}{x - x_0}.$$ Write a function that computes the derivative using the limit definition, using `limit`. # In[ ]: def lim_deriv(expr, x, x0): """ Computes the derivative of expr with respect to x at x0 using the limit definition. >>> lim_deriv(x**2, x, 0) 0 >>> lim_deriv(cos(x*y), x, pi) -y*sin(pi*y) Note that we must use this trick to take the derivative without evaluating at a point. >>> lim_deriv(exp(x**2), x, y).subs(y, x) 2*x*exp(x**2) """ # In[ ]: lim_deriv(x**2, x, 0) # In[ ]: lim_deriv(cos(x*y), x, pi) # The function you wrote above to compute limits using l'Hopital's rule is very fragile. And even if you try to make it sophisticated, it will still be unable to compute many limits. Try it on the following limits, and see what happens. Then try computing the same limits with `limit`. # # 1. $$\lim_{x\to 0}\frac{\log(x)}{x}$$ # 2. $$\lim_{x\to \infty}\frac{2^x}{3^x} \textbf{Warning: Be sure to save the notebook before you test this one, and be prepared to kill the kernel!}$$ # 3. $$\lim_{x\to \infty}x\sin{\left(\frac{1}{x}\right)}$$ # 4. $$\lim_{x\to 1}\arctan\left(\frac{1}{1 - x}\right)\; \text{Remember that $\arctan$ is called }\mathtt{atan}\text{ in SymPy}$$ # In[ ]: lhopital(log(x)/x, x, 0) # In[ ]: limit(log(x)/x, x, 0) # In[ ]: lhopital(2**x/3**x, x, oo) XXX: Don't run. This hangs the notebook # In[ ]: limit(2**x/3**x, x, oo) # In[ ]: lhopital(x**(1/x**2), x, 0) # In[ ]: limit(x**(1/x**.5), x, 0) # In[ ]: lhopital(x*sin(1/x), x, oo) # In[ ]: limit(x*sin(1/x), x, oo) # In[ ]: lhopital(atan(1/(1 - x)), x, 1) # In[ ]: limit(atan(1/(1 - x)), x, 1) # ## Series # The Fibonicci sequence is rexcursively defined by # # $$F_0 = 0,$$ # $$F_1 = 1,$$ # $$F_n = F_{n - 1} + F_{n - 2}.$$ # # The first few vales are 0, 1, 1, 2, 3, 5, 8, 13, 21, … # # The Fibonicci sequence has a generating function given by $$s(x) = \frac{x}{1 - x - x^2}$$ (see http://en.wikipedia.org/wiki/Fibonacci_number#Power_series for a derivation). What this means is that if we expand $s(x)$ as a power series, the coefficients are the Fibonicci numbers, $$s(x) = \sum_{n=0}^\infty F_nx^n$$ # Write a function that uses series to compute the nth Fibonicci number. # # Hint: `expr.coeff(x, n)` will give the coefficient of $x^n$ in an expression. For example # In[ ]: (1 + 2*x - x**2).coeff(x, 0) # In[ ]: (1 + 2*x - x**2).coeff(x, 1) # In[ ]: (1 + 2*x - x**2).coeff(x, 2) # In[ ]: def fib(n): """ Uses series expansion and a generating function to compute the nth Fibonnicci number. >>> fib(0) 0 >>> fib(4) 3 >>> fib(9) 34 """ # In[ ]: fib(0) # In[ ]: fib(4) # In[ ]: fib(9) # Note: if you really want to compute Fibonicci numbers, there is a function in SymPy called `fibonicci` that can do this far more efficiently. # In[ ]: [fibonacci(i) for i in range(10)] # `series` is nice if you want a fixed power series, but what if you don't know ahead of time how many terms you want? For that, there is the `lseries` method, which returns a generator of series terms. This is more efficient than recomputing the whole series again if you determine you need more terms. Here is an example usage (**Warning**: since series are in general infinite, `lseries` will return an infinite generator. Here we use `zip` to limit the number of terms). # In[ ]: for term, _ in zip(cos(x).lseries(x, 0), range(10)): display(term) # Write a function that computes the number of terms of a series expansion of a given function are needed to compute the given value near the given point within the given accuracy. For example, in the expansion of cos(x) near $\pi$, suppode we wish to compute $\pi + 1$. Let us see if terms up to $O(x^6)$ are sufficient (remember that due to a limitation in SymPy, when computing the series away from 0 with `series`, we have to use a shift) # In[ ]: a = cos(x).series(x, pi, 5).removeO().subs(x, x - pi).evalf(subs={x: pi + 1}) a # In[ ]: b = cos(pi + 1).evalf() b # So the expansion is accurate up to two places after the decimal point, i.e., within `0.01`. # In[ ]: abs(a - b) < 0.01 # Note, with `lseries`, we do not need to worry about shifts # In[ ]: for term, _ in zip(cos(x).lseries(x, pi), range(10)): display(term) # Hint: to get the exponent from a term like `3*(x - 1)**5`, use `as_coeff_exponent`. # In[ ]: (3*(x - 1)**5).as_coeff_exponent(x - 1) # In[ ]: def series_accuracy(func, expansion_point, evaluation_point, accuracy, x): """ Returns n such that series terms up to and including (x - expansion_point)**n (i.e., O((x - expansion_point)**(n + 1)) are needed to compute func at evaluation_point within the given accuracy. >>> series_accuracy(cos(x), pi, pi + 1, 0.01, x) 4 >>> series_accuracy(exp(x), 1, 10, 1, x) 23 """ # In[ ]: series_accuracy(cos(x), pi, pi + 1, 0.01, x) # In[ ]: series_accuracy(exp(x), 1, 10, 1, x) # In[ ]: