Lab Sheet 8

from header import *

Exercise 1

Here we define certain polynomials $p_0,\dotsc,p_3$.

x = sp.symbols('x')
p = [
 sp.sqrt(sp.Rational(1/2)),
 sp.sqrt(sp.Rational(3/2)) * x,
 sp.sqrt(sp.Rational(5/8)) * (3*x**2 - 1),
 sp.sqrt(sp.Rational(7/8)) * (5*x**3 - 3*x)
]
for i in range(4):
    display(Latex(f"$p_{i}= {sp.latex(p[i])}$"))

$p_0= \frac{\sqrt{2}}{2}$

$p_1= \frac{\sqrt{6} x}{2}$

$p_2= \frac{\sqrt{10} \cdot \left(3 x^{2} - 1\right)}{4}$

$p_3= \frac{\sqrt{14} \cdot \left(5 x^{3} - 3 x\right)}{4}$

We next calculate all the integrals $\int_{-1}^1 p_ip_j\,dx$ for $0\leq i,j<4$, and arrange them into a matrix. This turns out to be the identity matrix, which means that $\int_{-1}^1 p_ip_j\,dx$ is zero unless $i=j$, in which case it is equal to one.

sp.Matrix([[sp.integrate(p[i] * p[j],(x,-1,1)) for j in range(4)] for i in range(4)])

$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$

We now want to find a function of the form $p_4=ax^4+bx^2+c$ which makes the above pattern continue.

a,b,c = sp.symbols('a b c')
p4 = a * x**4 + b * x**2 + c
eqs = [sp.integrate(p4 * p[i],(x,-1,1)) for i in range(4)]
display(eqs)
eqs.append(sp.integrate(p4**2, (x,-1,1)) - 1)
sol = sp.solve(eqs)
display(sol)
p4 = p4.subs(sol[1])
p.append(p4)
display(Latex(f"$p_{4}= {sp.latex(p4)}$"))
display(Latex("This matrix of integrals should be the identity matrix"))
sp.Matrix([[sp.integrate(p[i] * p[j],(x,-1,1)) for j in range(5)] for i in range(5)])

[sqrt(2)*a/5 + sqrt(2)*b/3 + sqrt(2)*c,
 0,
 4*sqrt(10)*a/35 + 2*sqrt(10)*b/15,
 0]

[{a: -105*sqrt(2)/16, b: 45*sqrt(2)/8, c: -9*sqrt(2)/16},
 {a: 105*sqrt(2)/16, b: -45*sqrt(2)/8, c: 9*sqrt(2)/16}]

$p_4= \frac{105 \sqrt{2} x^{4}}{16} - \frac{45 \sqrt{2} x^{2}}{8} + \frac{9 \sqrt{2}}{16}$

This matrix of integrals should be the identity matrix

$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 1\end{matrix}\right]$

We now define $$ q(n) = \frac{\sqrt{n+1/2}}{2^n(n!)} \,\frac{d^n}{dx^n}((x^2-1)^n) $$

q = lambda n : sp.sqrt(n + sp.Rational(1,2))*sp.diff((x*x-1)**n,x,n)/(2**n * sp.factorial(n))

The cell below checks that $p_n=q(n)$ for $n\leq 4$. (This is Rodrigues's formula).

[sp.simplify(p[i] - q(i)) for i in range(5)]

[0, 0, 0, 0, 0]

del x, p, q

Exercise 2

Here we investigate integrals of the form $J(n,m) = \int x^n\ln(x)^m\,dx$.

x = sp.symbols('x')
def J(n,m):
    return sp.factor(sp.integrate(x **n * sp.log(x)**m, x))

Here are some examples:

display(Latex("$J(3,2) = " + sp.latex(J(3,2)) + "$"))
display(Latex("$J(4,1) = " + sp.latex(J(4,1)) + "$"))
display(Latex("$J(4,3) = " + sp.latex(J(4,3)) + "$"))

$J(3,2) = \frac{x^{4} \cdot \left(8 \log{\left(x \right)}^{2} - 4 \log{\left(x \right)} + 1\right)}{32}$

$J(4,1) = \frac{x^{5} \cdot \left(5 \log{\left(x \right)} - 1\right)}{25}$

$J(4,3) = \frac{x^{5} \cdot \left(125 \log{\left(x \right)}^{3} - 75 \log{\left(x \right)}^{2} + 30 \log{\left(x \right)} - 6\right)}{625}$

The pattern is that $J(n,m)$ is $x^{n+1}$ multiplied by a polynomial of degree $m$ in $\log(x)$.

del x

Exercise 3

x = sp.symbols('x')

Here are some integrals where sympy returns an answer in terms of functions that are not so familiar:

• $\text{erf}(x)$ is the error function
• $\text{li}(x)$ is the logarithmic integral function
• $F(x|a)$ is an elliptic integral
• $\Gamma(x)$ is the Gamma function
• ${}_2F_1(a,b,c|x)$ is a hypergeometric function

sp.integrate(sp.exp(-x*x), x)

$\displaystyle \frac{\sqrt{\pi} \operatorname{erf}{\left(x \right)}}{2}$

sp.integrate(1/sp.log(x), x)

$\displaystyle \operatorname{li}{\left(x \right)}$

sp.integrate(1/(sp.sqrt(1-x**2)*sp.sqrt(1-2*x**2)),x)

$\displaystyle \begin{cases} \frac{\sqrt{2} F\left(\operatorname{asin}{\left(\sqrt{2} x \right)}\middle| \frac{1}{2}\right)}{2} & \text{for}\: x > - \frac{\sqrt{2}}{2} \wedge x < \frac{\sqrt{2}}{2} \end{cases}$

sp.integrate((x**8+1)**sp.Rational(-1,2),x)

$\displaystyle \frac{x \Gamma\left(\frac{1}{8}\right) {{}_{2}F_{1}\left(\begin{matrix} \frac{1}{8}, \frac{1}{2} \\ \frac{9}{8} \end{matrix}\middle| {x^{8} e^{i \pi}} \right)}}{8 \Gamma\left(\frac{9}{8}\right)}$

Here are some examples that sympy cannot integrate, even in terms of obscure special functions like those mentioned above.

sp.integrate(sp.sin(x)*sp.log(sp.log(x)),x)

$\displaystyle \int \log{\left(\log{\left(x \right)} \right)} \sin{\left(x \right)}\, dx$

sp.integrate(sp.sin(sp.sin(sp.sin(x))),x)

$\displaystyle \int \sin{\left(\sin{\left(\sin{\left(x \right)} \right)} \right)}\, dx$

sp.integrate(1/sp.sqrt(1+x+x**10),x)

$\displaystyle \int \frac{1}{\sqrt{x^{10} + x + 1}}\, dx$

The simplest function that sympy cannot integrate is probably $x^x$:

display(sp.integrate(x+3,x))
display(sp.integrate(x*x,x))
display(sp.integrate(2**x,x))
display(sp.integrate(x**x,x))

$\displaystyle \frac{x^{2}}{2} + 3 x$

$\displaystyle \frac{x^{3}}{3}$

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

$\displaystyle \int x^{x}\, dx$

Exercise 4

Here we investigate the functions $y=\sum_{i=0}^3\sin((2i+1)x)/(2i+1)$ and $z=\int y\,dx$.

y = sp.sin(x) + sp.sin(3*x)/3 + sp.sin(5*x)/5 + sp.sin(7*x)/7
z = sp.integrate(y,x)
display(Latex("$y=" + sp.latex(y) + "$"))
display(Latex("$z=\\int y\\,dx=" + sp.latex(z) + "$"))

$y=\sin{\left(x \right)} + \frac{\sin{\left(3 x \right)}}{3} + \frac{\sin{\left(5 x \right)}}{5} + \frac{\sin{\left(7 x \right)}}{7}$

$z=\int y\,dx=- \cos{\left(x \right)} - \frac{\cos{\left(3 x \right)}}{9} - \frac{\cos{\left(5 x \right)}}{25} - \frac{\cos{\left(7 x \right)}}{49}$

We find that $y$ gives an approximate square wave, and $z$ gives an approximate sawtooth wave. (These are basic examples from Fourier analysis. The approximations would be better if we added more terms following the same pattern.)

xs = np.linspace(-3*np.pi, 3*np.pi,1000)
ys = sp.lambdify(x,y,'numpy')(xs)
zs = sp.lambdify(x,z,'numpy')(xs)
fig = go.Figure()
use_pi_ticks(fig,-3,3,1)
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='y'))
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='z'))

del x, y, z

Exercise 5

Here we investigate the functions $y=x\,e^{-x}\sin(20x)$ and $z=20\int y\,dx$.

x = sp.symbols('x')

y = x*sp.sin(20*x)*sp.exp(-x)
z = 20 * sp.integrate(y,x)
display(Latex("$y=" + sp.latex(y) + "$"))
display(Latex("$z=20\\int y\\,dx=" + sp.latex(z) + "$")) 
display(Latex("$\\qquad \\approx " + sp.latex(z.evalf()) + "$"))

$y=x e^{- x} \sin{\left(20 x \right)}$

$z=20\int y\,dx=- \frac{20 x e^{- x} \sin{\left(20 x \right)}}{401} - \frac{400 x e^{- x} \cos{\left(20 x \right)}}{401} + \frac{7980 e^{- x} \sin{\left(20 x \right)}}{160801} - \frac{800 e^{- x} \cos{\left(20 x \right)}}{160801}$

$\qquad \approx - 0.0498753117206983 x e^{- x} \sin{\left(20 x \right)} - 0.997506234413965 x e^{- x} \cos{\left(20 x \right)} + 0.0496265570487746 e^{- x} \sin{\left(20 x \right)} - 0.00497509343847364 e^{- x} \cos{\left(20 x \right)}$

The plots below show that the graph of $z$ is close to that of $y$, but shifted forwards by a quarter of a cycle. This is easy to understand in terms of the above approximate formula for $z$. Three of the numerical coefficients are small, and the remaining coefficient is close to $-1$, giving $$ z\approx -x\,e^{-x}\cos(20x)=x\,e^{-x}\sin(20x-\pi/2). $$

y_fun = sp.lambdify(x,y,'numpy')
z_fun = sp.lambdify(x,z,'numpy')
xs = np.linspace(0,10,1000)
ys = y_fun(xs)
zs = z_fun(xs)

fig = go.Figure()
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='y'))
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='z'))
move_legend(fig,0.5,0.9)
fig.show()

fig = make_subplots(rows=1, cols=2)
xs = np.linspace(0,3,1000)
ys = y_fun(xs)
zs = z_fun(xs)
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='y', line_color='blue'), row=1, col=1)
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='z', line_color='red'), row=1, col=1)
xs = np.linspace(1.255,1.57,1000)
ys = y_fun(xs)
zs = z_fun(xs)
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='y', line_color='blue', showlegend=False), row=1, col=2)
fig.add_trace(go.Scatter(x=xs, y=zs, mode='lines', name='z', line_color='red', showlegend=False), row=1, col=2)

We now consider $y^2+z^2$. In the approximate formula below, the terms $x^2e^{-2x}\sin^2(20x)$ and $x^2e^{-2x}\cos^2(20x)$ both appear with coefficients close to one, while all other coefficients are small. As $\sin^2(20x)+\cos^2(20x)=1$, this gives $y^2+z^2\approx x^2e^{-2x}$.

r2 = sp.simplify(y**2 + z**2)
r2_fun = sp.lambdify(x,r2,'numpy')
r2_approx = lambda x : x**2 * np.exp(-2*x)
display(Latex("$y^2+z^2=" + sp.latex(r2) + "$"))
display(Latex("$\\qquad\\approx " + sp.latex(sp.expand(r2.evalf())) + "$"))
xs = np.linspace(0,6,1000)
r2s = r2_fun(xs)
fig = go.Figure()
move_legend(fig,0.5,0.9)
fig.add_trace(go.Scatter(x=xs, y=r2s, mode='lines', name='$y^2+z^2$'))
fig.add_trace(go.Scatter(x=xs, y=r2_approx(xs), mode='lines', name='$x^2 e^{-2x}$'))

$y^2+z^2=\left(x^{2} \sin^{2}{\left(20 x \right)} + \frac{400 \left(401 x \sin{\left(20 x \right)} + 8020 x \cos{\left(20 x \right)} - 399 \sin{\left(20 x \right)} + 40 \cos{\left(20 x \right)}\right)^{2}}{25856961601}\right) e^{- 2 x}$

$\qquad\approx 1.00248754671924 x^{2} e^{- 2 x} \sin^{2}{\left(20 x \right)} + 0.0995018687694728 x^{2} e^{- 2 x} \sin{\left(20 x \right)} \cos{\left(20 x \right)} + 0.995018687694728 x^{2} e^{- 2 x} \cos^{2}{\left(20 x \right)} - 0.0049502800048653 x e^{- 2 x} \sin^{2}{\left(20 x \right)} - 0.0985093314251389 x e^{- 2 x} \sin{\left(20 x \right)} \cos{\left(20 x \right)} + 0.00992537344333894 x e^{- 2 x} \cos^{2}{\left(20 x \right)} + 0.00246279516451528 e^{- 2 x} \sin^{2}{\left(20 x \right)} - 0.000493793516694792 e^{- 2 x} \sin{\left(20 x \right)} \cos{\left(20 x \right)} + 2.47515547215435 \cdot 10^{-5} e^{- 2 x} \cos^{2}{\left(20 x \right)}$

del x, y, z, r2

Exercise 6

Here we find a function of the form $y=(ax^2+bx+c)e^{-x}$ such that $\int_{x=0}^\infty x^ky\,dx=0$ for $x=1,2,3$.

x, a, b, c = sp.symbols('x a b c')
f = (a*x**2+b*x+c)*sp.exp(-x)
sol = sp.solve([sp.integrate(x**k*f,(x,0,sp.oo))-k for k in range(1,4)])
f = f.subs(sol)
f

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

del x, a, b, c

Exercise 7

Here we find an approximate solution to $x(\ln(\ln(x))+\ln(x)-1)=541$ close to $x=100$.

x = sp.symbols('x')

y = x*(sp.log(sp.log(x)) + sp.log(x) - 1)
err = sp.lambdify(x,y-541,'numpy')

sc.optimize.root_scalar(err, bracket=[100,110]).root

104.35556573845588

del x

Exercise 8

Here we just want to generate the list $[a+b^2,a+b^3,a+b^4,\dotsc,a+b^{10}]$

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

[a + b**n for n in range(2, 11)]

[a + b**2,
 a + b**3,
 a + b**4,
 a + b**5,
 a + b**6,
 a + b**7,
 a + b**8,
 a + b**9,
 a + b**10]

del a, b

Exercise 9

Here we plot the curve $x^4+y^4=1$ together with the curve given parametrically by $x=\cos(t)/(2+\cos(4t))$ and $y=\sin(t)/(2+\cos(4t))$.

R=1.1
xs, ys = np.meshgrid(np.linspace(-R,R,100), np.linspace(-R,R,100))
zs = xs ** 4 + ys ** 4
ts = np.linspace(0,2*np.pi,200)
fig = go.Figure()
fig.update_layout(width=600, height=600, showlegend=False)
equal_aspect(fig)
fig.add_trace(go.Contour(x=xs[0], y=ys[:,0], z=zs, contours=dict(type='constraint', value=1), line_color='purple'))
fig.add_trace(go.Scatter(x=np.cos(ts)/(2 + np.cos(4*ts)), y=np.sin(ts)/(2 + np.cos(4*ts)), mode='lines', line_color='red'))

Exercise 10

Here we calculate $\pi^{76}/e^{87}$ to $100$ decimal places (showing that the result is quite close to $1$).

mpmath.mp.dps = 100
print(mpmath.pi**76/mpmath.exp(87))

0.9994714642786526084896787713200013244848205348814720991942990611898142176107601532266629087441688956

Exercise 11

Here we attempt to simplify the expression $y=\sin(2x)\tan(2x)\frac{d^2}{dx^2}\log(\tan(x))$.

x = sp.symbols('x')
y = sp.sin(2*x) * sp.tan(2*x) * sp.diff(sp.log(sp.tan(x)),x,2)
y

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

None of the following are very effective

display(sp.simplify(y))
display(sp.trigsimp(y))
display(sp.expand_trig(y))
display(sp.trigsimp(sp.expand_trig(y)))

$\displaystyle \frac{\left(\tan^{4}{\left(x \right)} - 1\right) \sin{\left(2 x \right)} \tan{\left(2 x \right)}}{\tan^{2}{\left(x \right)}}$

$\displaystyle \left(2 \tan^{2}{\left(x \right)} + 2 - \frac{8}{1 - \cos{\left(4 x \right)}}\right) \sin{\left(2 x \right)} \tan{\left(2 x \right)}$

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

$\displaystyle \frac{4 \cdot \left(\frac{2 \sin^{4}{\left(x \right)}}{\cos^{2}{\left(x \right)}} - 2 \tan^{2}{\left(x \right)} + 1\right)}{1 - 2 \cos^{2}{\left(x \right)}}$

However, sympy succeeds in showing that $y+4=0$, so we must have $y=-4$.

sp.simplify(y + 4)

$\displaystyle 0$

We could also test this by evaluating $y$ at some random values of $x$.

[y.evalf(subs={x:a}) for a in [-0.1, 0.4, 1.2, 3.1]]

[-4.00000000000000, -4.00000000000000, -4.00000000000000, -4.00000000000000]