Lab 1

Published

May 11, 2024

Goals for today

  1. Jump right in to learn enough of Python without any extra installation and without too much programming.
  2. Each student works through the exercise and picks up the skills needed to implement what you saw in the lecture in the computer.
  3. A self-evaluation form will be available about 15 minutes before the end of the class.

The start of every SymPy Live Shell

Go to SymPy Live Shell using your internet browser.

The live shell always runs the commands you see below at the beginning. This means that if you find yourself using SymPy elsewhere aside from the Live Shell, then you would expect some changes or even errors.

from sympy import *
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)

To appreciate the convenience of running these commands ahead without your intervention, try inputting the following commands one by one in the shell (and then pressing enter). Try to understand what happens.

a
x 
a + 1
x + 1

Our target is to use the shell interactively and gradually move on to using less and less of the shell. Once you exit the shell, everything you have done will be gone, ok?

But a very nice feature you should try is the “auto-complete” functionality, which is useful for not-so-short commands.

Learn to build symbolic expressions

  1. Operations: +, -, *, /, **
  2. Grouping symbols: ()

Some examples:

  1. \(x+y\) will appear if you type x+y in the shell.
  1. Try putting in some spaces, see what happens.
  2. Try producing an error by just typing x +, then pressing Enter.
  1. \(2xy\) is 2*x*y.
  1. Try inputting 2x*y instead. What happens and why do you think this happens?
  2. Suppose you want the symbolic expression \(2xa\). Try inputting 2*x*a. What happens and why do you think this happens?
  1. \(\dfrac{x}{y}\) is x/y.
  2. \(x^2\) is x**2. Take note that exponentiation is done differently than usual, so x^2 will NOT work.
  1. How about negative powers? What command should you enter to obtain \(x^{-2}\)? You may want to use grouping symbols.
  2. How about fractional powers? What command should you enter to obtain \(\sqrt[3]{x}\)? You may want to use grouping symbols.
  1. \(e^x\) is exp(x) and \(\ln x\) is log(x).
Note

Take note that \(\ln x=\log_e x\). But for some reason in economics, \(\ln x\) is sometimes written as \(\log x\) without the base \(e\).

  1. How about negative powers? What command should you enter to obtain \(x^{-2}\)? You may want to use grouping symbols.
  2. How about fractional powers? What command should you enter to obtain \(\sqrt[3]{x}\)? You may want to use grouping symbols.

You can actually produce the properties of logarithms in SymPy.

expand(log(x*y))

\(\displaystyle \log{\left(x y \right)}\)

What do you notice? Here we may have to incorporate assumptions about \(x\) and \(y\) to enable expansion.

x, y = symbols('x y', real = True, positive = True)
expand(log(x*y))

\(\displaystyle \log{\left(x \right)} + \log{\left(y \right)}\)

expand(log(x/y))

\(\displaystyle \log{\left(x \right)} - \log{\left(y \right)}\)

expand(log(x**2))

\(\displaystyle 2 \log{\left(x \right)}\)

What if you were asked to expand \(\ln(x^n)\)? How would you do it? Will it work if it were \(\ln(x^p)\)

Now, try practicing the input of expressions in the SymPy shell.

what would you input in the shell to produce

  1. \(\dfrac{t}{t^2+4}\) (Section 8.2, Example 1)
  2. \(e^{2x}-5e^x+4\) (Section 8.2, Example 2)
  3. \(e^{x-1}-x\) (Section 8.2, Example 3)
  4. \(\dfrac{\ln (3x^2+x^{1/3})}{2(3x+2)^3}\)

What if you want to input something like \(a^x\ln(x^b+c)\), for example? What do you have to do first?

What are the tasks we would ask the computer to do if we are in Section 8.2?

  1. Find stationary points. Solve equations in one variable.
  2. Determine which intervals would have a change in sign with respect to the first derivative.
  3. Make substitutions.
  4. Determine whether a function is convex or concave.

As long as we are clear about what tasks we want the shell to do, then it is really a matter of finding appropriate commands and syntax so that we would be able to carry out the tasks on the shell.

Consider Section 8.2, Example 1.

# Input the expression
expr = t/(t**2+4)
# Set up first derivative
Derivative(expr, t)

\(\displaystyle \frac{d}{d t} \frac{t}{t^{2} + 4}\)

# Find form of first derivative
Derivative(expr, t).doit()

\(\displaystyle - \frac{2 t^{2}}{\left(t^{2} + 4\right)^{2}} + \frac{1}{t^{2} + 4}\)

The preceding expression does not look exactly like what you see in the text. You have the option to simplify, expand, or factor. The portion _ refers to the most recently produced output. This saves a lot of typing.

_.expand()

\(\displaystyle - \frac{2 t^{2}}{t^{4} + 8 t^{2} + 16} + \frac{1}{t^{2} + 4}\)

_.simplify()

\(\displaystyle \frac{4 - t^{2}}{t^{4} + 8 t^{2} + 16}\)

_.factor()

\(\displaystyle - \frac{\left(t - 2\right) \left(t + 2\right)}{\left(t^{2} + 4\right)^{2}}\)

It seems directly using factor() works right away.

# Solve for stationary points
solve(Derivative(expr, t).doit(), t)
[-2, 2]

What you see here is a list in Python and the grouping symbols are square brackets. What this means is that \(t=-2\) and \(t=2\) solve \(c^{\prime}(t)=0\). We can rule out \(t=-2\) because of the condition in the problem where \(t\geq 0\).

To determine where \(c^\prime(t)\geq 0\) or \(c^\prime(t)\leq 0\), we have to look into how to solve inequalities.

# When is c'(t) greater than or equal to 0?
reduce_inequalities([Derivative(expr, t).doit() >= 0, t >= 0], t)

\(\displaystyle -2 \leq t \wedge 0 \leq t \wedge t \leq 2 \wedge t < \infty\)

The caret symbol ^ in the previous context has the meaning of AND in logic. Maybe it would be nicer to express the result in set notation:

reduce_inequalities([Derivative(expr, t).doit() >= 0, t >= 0], t).as_set()

\(\displaystyle \left[0, 2\right]\)

Pay attention to the formatting! Earlier you saw [-2, 2] and now you see \([0,2]\). Can you have an educated guess as to the difference?

Therefore, \(c^{\prime}(t)\geq 0\) whenever \(0\leq t\leq 2\).

Find when \(c^{\prime}(t)\leq 0\) by modifying the earlier commands.

Eventually, you will discover or at least be able to confirm the result found in the book. \(t=2\) maximizes \(c(t)\). The maximum value is \(c(2)\). To find this, we need to do a substitution:

# Extract the first solution in the list
t1 = solve(Derivative(expr, t).doit(), t)[0]
# Extract the second solution in the list
t2 = solve(Derivative(expr, t).doit(), t)[1]
# Symbolic result of substitution
expr.subs({t:t2})

\(\displaystyle \frac{1}{4}\)

# Numerical version
expr.subs({t:t2}).evalf()

\(\displaystyle 0.25\)

Notice that there is a new set of grouping symbols {}. This is called a dictionary in Python.

Although this was not discussed in Section 8.2 Example 1, you can choose to apply Theorem 8.2.2. We already know a stationary point. We just have to check whether \(c(t)\) is concave or convex over some interval \(I\).

# Set up second derivative
Derivative(expr, (t, 2))

\(\displaystyle \frac{d^{2}}{d t^{2}} \frac{t}{t^{2} + 4}\)

# Compute the second derivative
Derivative(expr, (t, 2)).doit()

\(\displaystyle \frac{2 t \left(\frac{4 t^{2}}{t^{2} + 4} - 3\right)}{\left(t^{2} + 4\right)^{2}}\)

Determine if it makes sense to simplify, expand, or factor the expression for the second derivative.

Now, is the second derivative positive or negative and over which intervals?

reduce_inequalities([_ >= 0, t >= 0], t).as_set()

\(\displaystyle \left\{0\right\} \cup \left[2 \sqrt{3}, \infty\right)\)

reduce_inequalities([Derivative(expr, (t, 2)).doit() <= 0, t >= 0], t).as_set()

\(\displaystyle \left[0, 2 \sqrt{3}\right]\)

What would be your conclusion?

Rough plots

plot(expr, title=None, xlabel=None, ylabel=None)

You can also focus on just \(t\geq 0\).

plot(expr, title=None, xlabel=None, ylabel=None, xlim = (0, 5))

You can choose to understand more of the behavior of the function so that you can improve the plot. For example, what happens to the function when \(t\to \pm\infty\), i.e., \[\lim_{t\to\infty} c(t)\,\, , \lim_{t\to -\infty} c(t)\]

Pay attention to how \(\infty\) is coded: it is two o’s or oo.

# Set up the one-sided limit as t goes to infinity
Limit(expr, t, oo)

\(\displaystyle \lim_{t \to \infty}\left(\frac{t}{t^{2} + 4}\right)\)

# Set up the one-sided limit as t goes to infinity
Limit(expr, t, oo).doit()

\(\displaystyle 0\)

Calculate the other limit for yourself.

Section 8.2 Exercise 10

Read the details of the exercise. The given \(f(x)\) is slightly different from the functions you have seen, like \(c(t)\) earlier. The function \(f(x)\) depends on parameters or unknown constants. Typically, there would be restrictions on these parameters and we could incorporate them as assumptions in symbolic computations.

# Set up the parameters of f(x)
I0, k, A, alpha = symbols('I_0 k A alpha', real = True, positive = True)

Try typing in I0 into the shell and press enter. What do you notice? Try it for alpha.

# Set up f(x)
expr = I0 + k*x + A*exp(-alpha*x)

Now, try to find minimizers of \(f(x)\). The result will then give you an answer to Section 8.2 Exercise 10a. You can test if the answers are indeed minimizers. After that, try (b).

Reporting results

Sometimes you may want the feeling that the shell is “talking” back to you. You ask it a question and then it gives you a response. Of course, you have to program the shell to do just that. The shell is not a living being (yet). But this would be a topic for next time. What is important now is for you to be able to use SymPy to implement what you see in ECOMATH.