Lab Sheet 5

from header import *

Exercise 1

(a)

The following lines define a "mystery function" $q(x)$.

qs = 'odpegd#{=#vs1ixqfwlrqv1vshfldo1huurubixqfwlrqv1ol+{,'
q = eval(''.join(list(map(lambda x: chr(ord(x)-3), qs))))

The function is undefined for $x<0$. It zero (and the graph is flat) at $x=0$, then it drops down to $-\infty$ at $x=1$, then climbs back to zero at about $x=1.4$, and increases thereafter towards $\infty$.

x = np.linspace(0,10,500)
y = [float(q(x0)) for x0 in x]
fig = go.Figure()
fig.update_layout(width=800, height=400)
fig.add_trace(go.Scatter(x=x, y=y, mode='lines'))

(b) When $x<0$, the function $q'(x)$ is undefined (because $q(x)$ is). We have $q'(0)=0$, and then $q'(x)<0$ for $0<x<1$ (because the graph is sloping downwards). The derivative is undefined again when $x=1$, but $q'(x)>0$ for $x>1$ (because the graph is sloping upwards).

(c)

def Q(x, h):
    return (q(x+h) - q(x))/h

print(Q(np.exp(2),0.01))
print(Q(np.exp(2),1e-5))
print(Q(np.exp(2),1e-10))

0.499830983358596
0.499999830783082
0.500000041370185

print(Q(np.exp(3),0.01))
print(Q(np.exp(3),1e-5))
print(Q(np.exp(3),1e-10))

0.333305681498786
0.333333305668759
0.333333360913457

print(Q(np.exp(4),0.01))
print(Q(np.exp(4),1e-5))
print(Q(np.exp(4),1e-10))

0.249994276886767
0.249999994395011
0.250004461577191

It is clear from the calculations above that $q'(e^{2})=\frac{1}{2}$ and $q'(e^{3})=\frac{1}{3}$ and $q'(e^{4})=\frac{1}{4}$. We therefore guess that $q'(e^{t})=\frac{1}{t}$ for all $t$. Moreover, given $x>0$, we can write $x=e^{t}$ with $t=\log(x)$, so $q'(x)=\frac{1}{\log(x)}$. This means that $q(x)$ is actually the same as a standard function called the logarithmic integral or $\text{li}(x)$. This is defined by $$ \text{li}(x)=\int_{t=0}^{x}\frac{dt}{\log(t)} . $$

del x
x, t = sp.symbols('x t')
sp.integrate(1/sp.log(t),(t,0,x)) - q(x)

$\displaystyle 0$

del x,t

Exercise 2

(a) The roots of $f'(x)$ occur where the graph of $f(x)$ is flat, or in other words, the tangent line is horizontal. This occurs wherever $f(x)$ has a local maximum or a local minimum (and possibly also in some other places, called inflexion points).

x = sp.symbols('x')

f = x**3 - x
display(Latex("$f(x) = " + sp.latex(f) + "$"))
df = sp.diff(f,x)
display(Latex("$f'(x) = " + sp.latex(df) + "$"))
f1 = sp.lambdify(x, f, 'numpy')
df1 = sp.lambdify(x, df, 'numpy')

$f(x) = x^{3} - x$

$f'(x) = 3 x^{2} - 1$

xs = np.linspace(-1.1,1.1,100)
ys = f1(xs)
roots = np.array([float(s[0]) for s in sp.solve([f],[x])])
fig = go.Figure()
fig.update_layout(width=600, height=400)
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$y=f(x)$'))
fig.add_trace(go.Scatter(x=roots, y=0*roots, mode='markers', marker=dict(color='red', size=8), name='roots'))
move_legend(fig, 0.7)
fig.show()
display(Latex(f"Roots of $f(x)$ are {roots}"))

Roots of $f(x)$ are [-1. 0. 1.]

(b)

xs = np.linspace(-1.1,1.1,100)
ys = f1(xs)
zs = df1(xs)
roots = np.array([float(s[0]) for s in sp.solve([df],[x])])
fig = go.Figure()
fig.update_layout(width=600, height=400, yaxis_range=[-1,1])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$f(x)$'))
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='$f\'(x)$'))
fig.add_trace(go.Scatter(x=roots, y=0*roots, mode='markers', marker=dict(color='red', size=8), showlegend=False))
fig.add_trace(go.Scatter(x=roots, y=f1(roots), mode='markers', marker=dict(color='blue', size=8), showlegend=False))
for r in roots:
    fig.add_trace(go.Scatter(x=[r,r], y=[0,f1(r)], mode='lines', line=dict(color='black', dash='dot'), showlegend=False))
move_legend(fig, 0.7)
fig.show()
display(Latex(f"Roots of $f'(x)$ are ${sp.latex(roots[0])}$ and ${sp.latex(roots[1])}$"))

Roots of $f'(x)$ are $-0.577350269189626$ and $0.577350269189626$

There are roots of $f(x)$ at $x=-1$, $x=0$ and $x=1$. Between $x=-1$ and $x=0$ we have a root of $f'(x)$ at about $x=-0.6$, corresponding to the top of the left-hand hump. Between $x=0$ and $x=1$ we have a root of $f'(x)$ at about $x=0.6$, corresponding to the top of the right-hand hump. Thus, Rolle's principle is satisfied.

(c)

g = sp.sin(x) + sp.sin(3*x)/3
display(Latex("$g(x) = " + sp.latex(g) + "$"))
dg = sp.diff(g,x)
display(Latex("$g'(x) = " + sp.latex(dg) + "$"))
g1 = sp.lambdify(x, g, 'numpy')
dg1 = sp.lambdify(x, dg, 'numpy')
xs = np.linspace(-2*np.pi,2*np.pi,400)
ys = g1(xs)
zs = dg1(xs)
fig = go.Figure()
fig.update_layout(width=800, height=400)
use_pi_ticks(fig, -2, 2)
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$g(x)$'))
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='$g\'(x)$'))

$g(x) = \sin{\left(x \right)} + \frac{\sin{\left(3 x \right)}}{3}$

$g'(x) = \cos{\left(x \right)} + \cos{\left(3 x \right)}$

We can see roots of $g(x)$ at $x=-2\pi$, $-\pi$, $0$, $\pi$ and $2\pi$; in general, the roots are at $x=n\pi$ for integers $n$. Between $x=0$ and $x=\pi$ there are two maxima and one local minimum, giving three roots of $g'(x)$. (This is perfectly consistent with Rolle's principle, which says only that there is at least one root of $g'(x) $ between $0$ and $\pi$.) Similarly, between $x=\pi$ and $x=2\pi$ there are two minima and one local maximum for $g(x)$, corresponding to the three places where the graph of $g'(x)$ (in orange) crosses the $x$-axis. We can find the location of the these roots as follows:

sp.solveset(dg)

$\displaystyle \left\{2 n \pi + \frac{\pi}{2}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{3 \pi}{2}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{5 \pi}{4}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{3 \pi}{4}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{7 \pi}{4}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{4}\; \middle|\; n \in \mathbb{Z}\right\}$

Sympy does not give the answer in the most efficient form: it is better to say that the roots are $(n+k/4)\pi$ with $n\in\mathbb{Z}$ and $k\in\{1,2,3\}$.

(d) Suppose that $a$ and $b$ are two different roots of $f(x)$, so $f(a)=0$ and $f(b)=0$. As $x$ moves from $a$ to $b$, the function $f(x)$ cannot increase all the time (otherwise $f(b)$ would be greater than 0) and it cannot decrease all the time (otherwise $f(b)$ would be less than 0). It must increase some of the time and decrease some of the time, so there must be some point at which it changes over from increasing to decreasing (or vice-versa), and at that point we will have $f'(x)=0$.

(e) As we see in the plot below, the roots of $h(x)=\tan(x)$ are at $x=n\pi$ for integers $n$.

h = sp.tan(x)
dh = sp.diff(h,x)
display(Latex(f"The derivative of $h(x) = {sp.latex(h)}$ is $h'(x) = {sp.latex(dh)}$"))
h1 = sp.lambdify(x, h, 'numpy')
dh1 = sp.lambdify(x, dh, 'numpy')
xs = np.linspace(-np.pi,np.pi,1000)
ys = h1(xs)
ix = np.abs(ys) < 10
xs = xs[ix]
ys = ys[ix]
ix0 = xs < -np.pi/2
ix1 = (-np.pi/2 < xs ) & (xs < np.pi/2)
ix2 = np.pi/2 < xs
fig = go.Figure()
fig.update_layout(width=800, height=400, yaxis_range=[-10,10], margin=dict(b=30))
use_pi_ticks(fig, -1, 1, 2)
# set figure title
fig.add_trace(go.Scatter(x=xs[ix0], y=ys[ix0], mode='lines', line=dict(color='blue'), name=r'$\tan(x)$'))
fig.add_trace(go.Scatter(x=xs[ix1], y=ys[ix1], mode='lines', line=dict(color='blue'), showlegend=False))
fig.add_trace(go.Scatter(x=xs[ix2], y=ys[ix2], mode='lines', line=dict(color='blue'), showlegend=False))
fig.add_trace(go.Scatter(x=[-np.pi/2, -np.pi/2], y=[-10, 10], mode='lines', line=dict(color='red', dash='dot'), showlegend=False))
fig.add_trace(go.Scatter(x=[ np.pi/2,  np.pi/2], y=[-10, 10], mode='lines', line=dict(color='red', dash='dot'), showlegend=False))

The derivative of $h(x) = \tan{\left(x \right)}$ is $h'(x) = \tan^{2}{\left(x \right)} + 1$

The graph is always sloping upwards, so $h'(x)>0$ for all $x$. In fact, we have $h'(x)=\sec(x)^{2}=1+\tan(x)^2$, which shows that $h'(x)\geq 1$ for all $x$. In any case, $h'(x)$ has no roots at all, contradicting Rolle's principle. So what goes wrong in the argument that we outlined? As $x$ moves from $0$ to $\pi$, the function increases from $0$ up towards $\infty$, then jumps down discontinuously. When $x=\pi/2$, neither $h(x)$ nor $h'(x) $ is really well-defined.

xs = np.linspace(-np.pi,np.pi,1000)
ys = dh1(xs)
ix = np.abs(ys) < 10
xs = xs[ix]
ys = ys[ix]
ix0 = xs < -np.pi/2
ix1 = (-np.pi/2 < xs ) & (xs < np.pi/2)
ix2 = np.pi/2 < xs

fig = go.Figure()
fig.update_layout(width=800, height=400, yaxis_range=[0,10], margin=dict(t=50, b=30))
use_pi_ticks(fig, -1, 1, 2)
fig.update_layout(title=r"$\tan'(x) = 1 + \tan^2(x)$")
fig.add_trace(go.Scatter(x=xs[ix0], y=ys[ix0], mode='lines', line=dict(color='blue'), showlegend=False))
fig.add_trace(go.Scatter(x=xs[ix1], y=ys[ix1], mode='lines', line=dict(color='blue'), showlegend=False))
fig.add_trace(go.Scatter(x=xs[ix2], y=ys[ix2], mode='lines', line=dict(color='blue'), showlegend=False))
fig.add_trace(go.Scatter(x=[-np.pi/2, -np.pi/2], y=[0, 10], mode='lines', line=dict(color='red', dash='dot'), showlegend=False))
fig.add_trace(go.Scatter(x=[ np.pi/2,  np.pi/2], y=[0, 10], mode='lines', line=dict(color='red', dash='dot'), showlegend=False))

del f,g,h,df,dg,dh,f1,g1,h1,df1,dg1,dh1,xs,ys,zs,ix,ix0,ix1,ix2

Exercise 3

x, n, a, b, c, d, p = sp.symbols('x n a b c d p')

def schwarzian(f,x):
    return f.diff(x,3)/f.diff(x) - sp.Rational(3,2) * (f.diff(x,2)/f.diff(x))**2

Put $y=x^{n}$. Then $y'=nx^{n-1}$ and $y''=n(n-1) x^{n-2}$ and $ {y'''}=n (n-1) (n-2) x^{n-3} $. If $n=0$ then $y'=0$, so we cannot divide by $y'$, so $S(y)$ is undefined. For any other value of $n$, we have $y''/y'=(n-1)/x$ and $y'''/y'=(n-1)(n-2)/x^{2}$. This gives \begin{align*} S(y) &= \frac{{y'''}}{{y'}}-\frac{3 (\frac{{y''}}{{y'}})^{2}}{2} =\frac{(n -1) (n -2)}{x^{2}}-\frac{3 (\frac{n -1}{x})^{2}}{2} \\ &=\frac{2 (n^{2}-3 n +2)-3 (n -1)^{2}}{2 x^{2}} \\ &=\frac{1-n^{2}}{2 x^{2}} \end{align*}

This is zero for $n=1$ or $n=-1$. We can also do the calculation using sympy as follows:

ans = sp.simplify(schwarzian(x ** n, x))
display(ans)
display([s[0] for s in sp.solve([ans],[n])])

$\displaystyle \frac{1 - n^{2}}{2 x^{2}}$

[-1, 1]

y = (a * x + b) / (c * x + d)
y1 = sp.factor(sp.expand(y.diff(x)))
y2 = sp.factor(sp.expand(y1.diff(x)))
y3 = sp.factor(sp.expand(y2.diff(x)))
Sy = sp.simplify(schwarzian(y,x))

display(Latex(
  r"\begin{align*}" + 
  f"y &= {sp.latex(y)} \\\\" + 
  f"y' &= {sp.latex(y1)} \\\\" + 
  f"y'' &= {sp.latex(y2)} \\\\" + 
  f"y''' &= {sp.latex(y3)} \\\\" +
  f"S(y) &= {sp.latex(Sy)}" 
  r"\end{align*}"
))

\begin{align*}y &= \frac{a x + b}{c x + d} \\y' &= \frac{a d - b c}{\left(c x + d\right)^{2}} \\y'' &= - \frac{2 c \left(a d - b c\right)}{\left(c x + d\right)^{3}} \\y''' &= \frac{6 c^{2} \left(a d - b c\right)}{\left(c x + d\right)^{4}} \\S(y) &= 0\end{align*}

z = y.subs({x : p ** x})
Sz = sp.simplify(schwarzian(z,x))
Sz

$\displaystyle - \frac{\log{\left(p \right)}^{2}}{2}$

def T(y, x):
    return sp.sqrt(sp.diff(y,x)) * sp.diff(1/sp.sqrt(sp.diff(y,x)),x,2)

y = x ** n
display(sp.simplify(T(y, x))) 
display(sp.simplify(schwarzian(y,x)))

$\displaystyle \frac{n^{2} - 1}{4 x^{2}}$

$\displaystyle \frac{1 - n^{2}}{2 x^{2}}$

y = (a * x + b)/(c * x + d)
display(sp.simplify(T(y, x))) 
display(sp.simplify(schwarzian(y,x)))

$\displaystyle 0$

$\displaystyle 0$

y = sp.log(x)
display(sp.simplify(T(y, x))) 
display(sp.simplify(schwarzian(y,x)))

$\displaystyle - \frac{1}{4 x^{2}}$

$\displaystyle \frac{1}{2 x^{2}}$

u = sp.Function('u')
sp.simplify(2 * T(u(x),x) + schwarzian(u(x),x))

$\displaystyle 0$

del x, n, a, b, c, d, p, y, y1, y2, y3, Sy, z, Sz, T, u

Exercise 4

x = sp.symbols('x')
y = sp.cos(-2 * sp.log(x))
dydx = sp.diff(y, x)
display(Latex(f"$y = {sp.latex(y)}$"))
display(Latex(f"$y' = {sp.latex(dydx)}$"))
f = sp.lambdify(x, y, 'numpy')
f1 = sp.lambdify(x, dydx, 'numpy')

$y = \cos{\left(2 \log{\left(x \right)} \right)}$

$y' = - \frac{2 \sin{\left(2 \log{\left(x \right)} \right)}}{x}$

fig = make_subplots(rows=2, cols=1)
fig.update_layout(width=800, height=800)
xs = np.linspace(1e-6,0.02,1000)
fig.add_trace(go.Scatter(x=xs, y=f(xs), mode='lines', name=r'$y=\cos(-2\log(x))$'), row=1, col=1)
xs = np.linspace(1e-6,5000,1000)
fig.add_trace(go.Scatter(x=xs, y=f(xs), mode='lines', name=r'$y=\cos(-2\log(x))$'), row=2, col=1)

ks = np.arange(10)
a = np.exp((1/4 - ks/2) * np.pi)
m = f1(a)

k = 9
fig = go.Figure()
fig.update_layout(width=1000, height=400)
xs = np.linspace(10**(-k-3),2*a[k],1000)
fig.add_trace(go.Scatter(x=xs, y=f(xs), mode='lines', showlegend=False))
fig.add_trace(go.Scatter(x=xs, y=m[k] * (xs - a[k]), mode='lines', line=dict(color='orange'), showlegend=False))
fig.add_trace(go.Scatter(x=[a[k]],y=[0], mode='markers', marker=dict(color='blue', size=8), showlegend=False))
fig.add_trace(go.Scatter(x=[0],y=[2 * (-1)**k], mode='markers', marker=dict(color='red', size=8), showlegend=False))