Lab Sheet 6

from header import *

Exercise 1.1

Here are the graphs of $y=x^{n}e^{-x}$ for $n =3,4$ and $5$:

xs = np.linspace(0, 10, 200)
fig = make_subplots(rows=1, cols=3)
fig.update_layout(width=1200, height=400)
for i, n in enumerate(range(3, 6)):
    fig.add_trace(
        go.Scatter(x=xs, y=xs ** n * np.exp(-xs), mode='lines', showlegend=False),
        row=1, col=i + 1
    )
    fig.update_xaxes(dtick=1, row=1, col=i + 1)
    fig.add_annotation(text=f"$x^{{{n}}} e^{{-x}}$",x=8, y = 0.8, yref=f'y{n+1} domain', row=1, col=i + 1,showarrow=False, font=dict(size=20))

fig.show()

You can just about see that the peaks occur at $x=3$, $x=4$ and $x=5$, suggesting that for general $n$ we should have a peak at $x=n$. To check this, we must solve $dy/dx=0$.

x = sp.symbols('x')
n = sp.symbols('n', integer=True, positive=True)
y = x ** n * sp.exp(-x)
y1 = sp.simplify(sp.diff(y, x))
x0 = sp.solve([y1],[x])[1][0]
display(Latex(f"The derivative of ${sp.latex(y)}$ is ${sp.latex(y1)}$, and this vanishes at $x={sp.latex(x0)}$"))

The derivative of $x^{n} e^{- x}$ is $x^{n - 1} \left(n - x\right) e^{- x}$, and this vanishes at $x=n$

y_max = y.subs(x, x0)
display(Latex("The maximum value of $" + sp.latex(y) + "$ is $y_{\\text{max}}=" + sp.latex(y_max) + "$"))
z = sp.simplify(y / y_max)
z_fun = sp.lambdify((n, x), z, 'numpy')
display(Latex("The normalised function $z=y/y_{\\text{max}}$ is given by $" + sp.latex(z) + "$"))

The maximum value of $x^{n} e^{- x}$ is $y_{\text{max}}=n^{n} e^{- n}$

The normalised function $z=y/y_{\text{max}}$ is given by $\left(\frac{x}{n}\right)^{n} e^{n - x}$

We can plot $z$ for different values of $n$ as shown below.

xs = np.linspace(0,20,500)
fig = go.Figure()
for n0 in range(1,9):
    fig.add_trace(go.Scatter(x=xs, y=z_fun(n0, xs), mode='lines', name=f'$n={n0}$'))

fig.show()

del x

Exercise 1.2

x = sp.symbols('x')
p = -10*x**6+156*x**5-945*x**4+2780*x**3-4080*x**2+2880*x
display(Latex("We consider the function $p(x)=" + sp.latex(p) + "$"))
p1 = sp.diff(p, x)
display(Latex("The derivative of is $p'(x)=" + sp.latex(p1) + "=" + sp.latex(sp.factor(p1)) + "$"))
p2 = sp.diff(p1, x)
p_fun = sp.lambdify(x, p, 'numpy')

root_multiplicities = sp.roots(sp.diff(p, x))
roots = sorted(list(root_multiplicities.keys()))
root_multiplicities = [root_multiplicities[r] for r in roots]
display(Latex(f"The roots of $p'(x)$ are {roots} with multiplicities {root_multiplicities}"))

We consider the function $p(x)=- 10 x^{6} + 156 x^{5} - 945 x^{4} + 2780 x^{3} - 4080 x^{2} + 2880 x$

The derivative of is $p'(x)=- 60 x^{5} + 780 x^{4} - 3780 x^{3} + 8340 x^{2} - 8160 x + 2880=- 60 \left(x - 4\right)^{2} \left(x - 3\right) \left(x - 1\right)^{2}$

The roots of $p'(x)$ are [1, 3, 4] with multiplicities [2, 1, 2]

xs = np.linspace(0,5.2023,1000)
fig = go.Figure()
fig.update_layout(width=800, height=400, showlegend=False)
fig.add_trace(go.Scatter(x=xs, y=p_fun(xs), mode='lines', name='$p(x)$'))
for a in roots:
    a0 = float(a)
    b0 = float(p_fun(a0))
    fig.add_trace(go.Scatter(x=[a0], y=[b0], mode='markers', name=f'$p({a})={b}$', marker=dict(color='red', size=10)))
    fig.add_trace(go.Scatter(x=[a0-0.3,a0+0.3], y=[b0,b0], mode='lines', line=dict(color='green')))
    fig.add_trace(go.Scatter(x=[a0,a0], y=[0,b0], mode='lines', line=dict(color='green', dash='dash')))

fig.show()

for a in roots:
    b = p.subs(x, a)
    c = p2.subs(x, a)
    if c < 0:
       t = 'a local maximum'
    elif c > 0:
       t = 'a local minimum'
    else:
       t = 'an inflexion point'
    display(Latex(f"At $x={a}$ we have $p(x)={b}$ and $p''(x)={c}$.  This is {t}."))

At $x=1$ we have $p(x)=781$ and $p''(x)=0$. This is an inflexion point.

At $x=3$ we have $p(x)=1053$ and $p''(x)=-240$. This is a local maximum.

At $x=4$ we have $p(x)=1024$ and $p''(x)=0$. This is an inflexion point.

del x, p, p1, p2, p_fun, a, b, c, roots, root_multiplicities

Exercise 2.1

x, y, t = sp.symbols('x y t')
u = x * sp.sin(x**2 + y**2) - y * sp.cos(x**2 + y**2)
u_fun = sp.lambdify((x, y), u, 'numpy')

display(Latex("We plot the curve given by $u(x,y)=0$, where $u=" + sp.latex(u) + "$."))
xs, ys = np.meshgrid(np.linspace(-5,5,200), np.linspace(-5,5,200))
zs = u_fun(xs, ys)
fig = go.Figure()
fig.update_layout(width=600, height=600, xaxis_range=[-4,4], yaxis_range=[-4,4], showlegend=False)
fig.add_trace(go.Contour(x=xs[0], y=ys[:,0], z=zs, contours=dict(type='constraint', value=0), line_color='purple'))

We plot the curve given by $u(x,y)=0$, where $u=x \sin{\left(x^{2} + y^{2} \right)} - y \cos{\left(x^{2} + y^{2} \right)}$.

The slope of the curve can be given in terms of $x$ and $y$ by implicit differentiation.

slope1 = sp.idiff(u, y, x)
slope1

$\displaystyle - \frac{2 x^{2} \cos{\left(x^{2} + y^{2} \right)} + 2 x y \sin{\left(x^{2} + y^{2} \right)} + \sin{\left(x^{2} + y^{2} \right)}}{2 x y \cos{\left(x^{2} + y^{2} \right)} + 2 y^{2} \sin{\left(x^{2} + y^{2} \right)} - \cos{\left(x^{2} + y^{2} \right)}}$

The above expression can be simplified using the fact that $u(x,y)=0$, which gives $\sin(x^2+y^2)=yx^{-1}\cos(x^2+y^2)$. After making this substitution, all the trigonometric functions cancel out.

slope1 = sp.simplify(slope1.subs(sp.sin(x**2+y**2), y/x*sp.cos(x**2+y**2)))
slope1

$\displaystyle \frac{- 2 x \left(x^{2} + y^{2}\right) - y}{x \left(2 x y - 1\right) + 2 y^{3}}$

xt = t * sp.cos(t**2)
yt = t * sp.sin(t**2)
display(Latex("The curve can also be described parametrically by $x(t)=" + sp.latex(xt) + "$ and $y(t)=" + sp.latex(yt) + "$."))
display(Latex("To verify this, we just need to check that $u(x(t),y(t))=0$."))
sp.simplify(u.subs({x : xt, y : yt}))

The curve can also be described parametrically by $x(t)=t \cos{\left(t^{2} \right)}$ and $y(t)=t \sin{\left(t^{2} \right)}$.

To verify this, we just need to check that $u(x(t),y(t))=0$.

$\displaystyle 0$

The plot below is based on the parametric description; it is visibly the same as the previous plot.

xy_fun = sp.lambdify(t, [xt,yt], 'numpy')
ts = np.linspace(-5,5,1000)
xys = xy_fun(ts)
fig = go.Figure()
fig.update_layout(width=600, height=600, xaxis_range=[-4,4], yaxis_range=[-4,4], showlegend=False)
fig.add_trace(go.Scatter(x=xys[0], y=xys[1], mode='lines', line_color='blue'))

From the parametric description we obtain another expression for the slope:

slope2 = yt.diff(t)/xt.diff(t)
slope2

$\displaystyle \frac{2 t^{2} \cos{\left(t^{2} \right)} + \sin{\left(t^{2} \right)}}{- 2 t^{2} \sin{\left(t^{2} \right)} + \cos{\left(t^{2} \right)}}$

The following calculation verifies that the two formulae for the slope are equivalent:

sp.simplify(slope1.subs({x : xt, y : yt}) - slope2)

$\displaystyle 0$

del x, y, t, u, u_fun, xs, ys, zs, xt, yt, xy_fun, ts, xys, slope1, slope2

Exercise 2.2

x, y, t = sp.symbols('x y t')
u = (x**2 + y**2) ** 2 + 85 * (x**2 + y**2) - 500 + 18*x*(3*y**2 - x**2)
u_fun = sp.lambdify((x, y), u, 'numpy')

xt = 6*sp.cos(t) + 8*sp.cos(t)**2 - 4
yt = 2*sp.sin(t)*(3-4*sp.cos(t))
xy_fun = sp.lambdify(t, [xt,yt], 'numpy')

display(Latex("On the left we plot the curve given by $u(x,y)=0$, where $u=" + sp.latex(u) + "$."))
display(Latex("On the right we plot the curve given parametrically by $x(t)=" + sp.latex(xt) + "$ and $y(t)=" + sp.latex(yt) + "$."))

xs, ys = np.meshgrid(np.linspace(-10,10,200), np.linspace(-10,10,200))
zs = u_fun(xs, ys)
xys = xy_fun(np.linspace(0,2*np.pi,1000))

fig = make_subplots(rows=1, cols=2, subplot_titles=['Implicit plot', 'Parametric plot'])
fig.update_layout(width=1000, height=500, showlegend=False)
for c in range(1,3):
    fig.update_xaxes(range=[-10,10], row=1, col=c)
    fig.update_yaxes(range=[-10,10], row=1, col=c)
fig.add_trace(go.Contour(x=xs[0], y=ys[:,0], z=zs, contours=dict(type='constraint', value=0), line_color='purple'),row=1, col=1)
fig.add_trace(go.Scatter(x=xys[0], y=xys[1], mode='lines', line_color='blue'), row=1, col=2)

On the left we plot the curve given by $u(x,y)=0$, where $u=85 x^{2} + 18 x \left(- x^{2} + 3 y^{2}\right) + 85 y^{2} + \left(x^{2} + y^{2}\right)^{2} - 500$.

On the right we plot the curve given parametrically by $x(t)=8 \cos^{2}{\left(t \right)} + 6 \cos{\left(t \right)} - 4$ and $y(t)=2 \left(3 - 4 \cos{\left(t \right)}\right) \sin{\left(t \right)}$.

The following calculation shows that $u(x(t),y(t))=0$, so the two curves are the same.

sp.simplify(sp.expand(u.subs({x : xt, y : yt})))

$\displaystyle 0$

slope1 = sp.simplify(sp.idiff(u, y, x))
display(Latex("Implicit formula for the slope: $$" + sp.latex(slope1) + "$$"))
slope2 = yt.diff(t)/xt.diff(t)
display(Latex("Parametric formula for the slope: $$" + sp.latex(slope2) + "$$"))
err = sp.trigsimp(sp.factor(slope1.subs({x : xt, y : yt}) - slope2))
display(Latex("Consistency check: $" + sp.latex(err) + "$"))

Implicit formula for the slope: $$\frac{- 2 x^{3} + 27 x^{2} - 85 x - y^{2} \cdot \left(2 x + 27\right)}{y \left(2 x^{2} + 54 x + 2 y^{2} + 85\right)}$$

Parametric formula for the slope: $$\frac{2 \cdot \left(3 - 4 \cos{\left(t \right)}\right) \cos{\left(t \right)} + 8 \sin^{2}{\left(t \right)}}{- 16 \sin{\left(t \right)} \cos{\left(t \right)} - 6 \sin{\left(t \right)}}$$

Consistency check: $0$

Exercise 3.1

We define $r(n) = \frac{d^n}{dx^n}\left(\frac{x^n\log(x)}{n!}\right)$

x = sp.symbols('x')

def r(n):
    return sp.diff(x ** n * sp.log(x) / sp.factorial(n),x,n) 

The list $r(0),\dotsc,r(10)$ is as follows:

rr = [r(n) for n in range(0,11)]
rr

[log(x),
 log(x) + 1,
 log(x) + 3/2,
 log(x) + 11/6,
 log(x) + 25/12,
 log(x) + 137/60,
 log(x) + 49/20,
 log(x) + 363/140,
 log(x) + 761/280,
 log(x) + 7129/2520,
 log(x) + 7381/2520]

It seems experimentally that $r(n)-r(n-1)=1/n$, so $r(n)=x+\sum_{k=1}^n\frac{1}{k}$.

rr_dif = [r(n) - r(n-1) for n in range(1,11)]
rr_dif

[1, 1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10]

The sum can be expressed in terms of

• A certain constant called sp.EulerGamma in scipy (or just $\gamma$ in traditional notation)
• A certain function called sp.digamma(x) in scipy (or $\Psi(x)$ in traditional notation)

def r1(n):
    return sp.log(x) + sp.digamma(n+1) + sp.EulerGamma

The calculation below shows that r(n)=r1(n) for $n=1,\dotsc,10$. In fact this holds for all $n$, as one can prove by induction, but we will not give the details here.

[r(n) - r1(n) for n in range(1,11)]

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

del x, r, rr, rr_dif, r1

Exercise 3.2

Here we aim to find a differential equation that is satisfied by the function $y=t^2 e^t$. More precisely, we want to find constants $a$, $b$ and $c$ such that $y'''+ay''+by'+cy=0$.

We first take $a$, $b$ and $c$ to be symbols, and calculate $y'''+ay''+by'+cy$.

t, a, b, c = sp.symbols('t a b c')

y = t ** 2 * sp.exp(t)
display(y)
z = sp.simplify(sp.expand(sp.diff(y,t,3) + a * sp.diff(y,t,2) + b * sp.diff(y,t) + c * y))
display(z)

$\displaystyle t^{2} e^{t}$

$\displaystyle \left(a t^{2} + 4 a t + 2 a + b t^{2} + 2 b t + c t^{2} + t^{2} + 6 t + 6\right) e^{t}$

The polynomial part of $z$ can be accessed as z.args[0]. We extract the list of coefficients of this polynomial and then use sp.solve() to find values of $a$, $b$ and $c$ which make these coefficients equal to zero.

sol = sp.solve(sp.Poly(z.args[0],t).coeffs())
sol

{a: -3, b: 3, c: -1}

The conclusion is that $y'''-3y''+3y'-y=0$. (The coefficients are the same as in the binomial expansion of $(u-1)^3$. This is not a coincidence, but we will not explain the reason here.)

del t, a, b, c, y, z, sol

Exercise 3.3

We define $$ p(n) = e^{x^2} \frac{d^n}{dx^n}e^{-x^2}. $$

x = sp.symbols('x')

def p(n):
    return sp.expand(sp.exp(x**2)*sp.diff(sp.exp(-x**2),x,n))

We calculate $p(n)$ for $n=0,\dotsc,10$.

for n in range(0,11):
    display(p(n))

$\displaystyle 1$

$\displaystyle - 2 x$

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

$\displaystyle - 8 x^{3} + 12 x$

$\displaystyle 16 x^{4} - 48 x^{2} + 12$

$\displaystyle - 32 x^{5} + 160 x^{3} - 120 x$

$\displaystyle 64 x^{6} - 480 x^{4} + 720 x^{2} - 120$

$\displaystyle - 128 x^{7} + 1344 x^{5} - 3360 x^{3} + 1680 x$

$\displaystyle 256 x^{8} - 3584 x^{6} + 13440 x^{4} - 13440 x^{2} + 1680$

$\displaystyle - 512 x^{9} + 9216 x^{7} - 48384 x^{5} + 80640 x^{3} - 30240 x$

$\displaystyle 1024 x^{10} - 23040 x^{8} + 161280 x^{6} - 403200 x^{4} + 302400 x^{2} - 30240$

It seems that $p(n)$ is always a polynomial in $x$ with highest term $(-2x)^n$.

Next, note that when $n$ is even we only get even powers of $x$, and when $n$ is odd we only get odd powers of $x$. For example, $p(7)$ involves only $x,x^{3},x^{5}$ and $x^{7}$, whereas $p (6)$ involves $x^{2},x^{4}$ and $x^{6}$ (and a constant term, which we can think of as a multiple of $x^{0}$).

Now look at the last term in $p(n)$. When $n$ is even, the last term is a constant, but when $n$ is odd, it is a multiple of $x$. It is best to consider these separately, starting with the even case, where we may write $n=2m$.

We see that the last term is $(-1)^{m}$ times a constant, with the sequence of constants being $1,2,12,120,1680,30240,665280$ and so on. If we enter this in the Online Encyclopedia of Integer Sequences we get an answer including the line

A001813: Quadruple factorial numbers: a(n) = (2n)!/n!

[p(2*m).subs({x : 0}) for m in range(10)]

[1, -2, 12, -120, 1680, -30240, 665280, -17297280, 518918400, -17643225600]

We can now check directly that the constant term in $p(2m)$ is $(-1)^m(2m)!/m!$, at least for $m<10$:

[(-1)**m * sp.factorial(2*m)/sp.factorial(m) for m in range(10)]

[1, -2, 12, -120, 1680, -30240, 665280, -17297280, 518918400, -17643225600]

Of course this kind of experimental approach does not really prove that the formula is correct, but it is very suggestive.