Lab Sheet 7

from header import *

Integrals of rational functions can be expressed in a number of different, but mathematically equivalent, ways. The standard function sp.integral(f, x) returns a result which depends in a rather random way on the detailed structure of the function f. We therefore write a function ratint(f, x) that expresses the integral in a more coherent and systematic way.

def ratint(f, x):
    f = sp.factor(f)
    nf = sp.expand(sp.numer(f))
    df = sp.expand(sp.denom(f))
    nfd = sp.degree(nf, x)
    dfd = sp.degree(df, x)
    d = nfd - dfd
    rs0 = sp.Poly(df,x).all_roots(multiple=False)
    rs = []
    for r in rs0:
        if isinstance(r[0],sp.CRootOf):
            r = (r[0].evalf(),r[1])
        rs.append(r)
    c = df.coeff(x, dfd)

    i = 0
    ps = []
    for r, m in rs:
        if sp.im(r) < 0:
            continue
        elif sp.im(r) == 0:
            ps.append((i, x - r, m))
            i += 1
        else:
            ps.append((i, x**2 - 2*sp.re(r)*x + sp.re(r)**2 + sp.im(r)**2, m))
            i += 1
    qs = []
    for i, p in enumerate(ps):
        qs.append((p[0], p[1], p[2], c*sp.Mul(*[p1[1] ** p1[2] for p1 in ps if p1[0] != p[0]])))
    ms = [[x**i, x**i*df, x**(i+1)/(i+1)] for i in range(d+1)]
    for q in qs:
        n = sp.degree(q[1], x)
        m = q[2]
        if n == 1:
            ms.append([1/q[1],q[1]**(m-1)*q[3],sp.log(abs(q[1]))])
            for i in range(2, m+1):
                ms.append([q[1]**(-i), q[1]**(m-i)*q[3], q[1]**(1-i)/(1-i)])
        else:
            a = q[1].coeff(x, 1)
            b = q[1].coeff(x, 0)
            u = 1/sp.sqrt(b - a**2/4)
            v = a*u/2
            ms.append([1/q[1], q[1]**(m-1)*q[3], u * sp.atan(u*x+v)])
            ms.append([(2*x+a)/q[1], (2*x+a) * q[1]**(m-1)*q[3], sp.log(abs(q[1]))])
            for i in range(1, 2*m-1):
                u = i * x**(i-1) * q[1] + (1-m) * x**i * (2*x+a)
                ms.append([u*q[1]**(-m), u, x**i*q[1]**(1-m)])

    a = sp.IndexedBase('a')
    err = sp.collect(sp.expand(nf - sp.Add(*[a[i] * m[1] for i, m in enumerate(ms)])),x)
    sol = sp.solve(sp.Poly(err,x).coeffs())
    F = sp.Add(*[a[i] * m[2] for i, m in enumerate(ms)]).subs(sol)
    return F

We now define a list of functions $g_0,\dotsc,g_9$ which are good test cases for our integration routine.

x = sp.symbols('x')

g = [
 x**2 + x + 1 + 1/x + 1/x**2,
 1/(x-1) + 1/(x**2-1) + 1/(x**3-1),
 1/(x-1)**2 + 1/(x-2)**3 + 1/(x-3)**4,
 1/(x+1) + 1/(x**2+1) + 1/(x**3+1),
 x/(x+1)*(x+2)/(x+3)*(x+4),
 x/(x+1)**2*(x+2)/(x+3)**2*(x+4),
 (x**3+4*x**2+x+2)/(x-1)**2/(x+1)**3,
 1/(x**2+8*x+1),
 (x**2+1)**(-5),
 1/(x**8+x**2+1),
]

for i, g0 in enumerate(g):
    display(Latex(f"$g_{i}(x) = {sp.latex(g0)}$"))

$g_0(x) = x^{2} + x + 1 + \frac{1}{x} + \frac{1}{x^{2}}$

$g_1(x) = \frac{1}{x^{3} - 1} + \frac{1}{x^{2} - 1} + \frac{1}{x - 1}$

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

$g_3(x) = \frac{1}{x^{3} + 1} + \frac{1}{x^{2} + 1} + \frac{1}{x + 1}$

$g_4(x) = \frac{x \left(x + 2\right) \left(x + 4\right)}{\left(x + 1\right) \left(x + 3\right)}$

$g_5(x) = \frac{x \left(x + 2\right) \left(x + 4\right)}{\left(x + 1\right)^{2} \left(x + 3\right)^{2}}$

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

$g_7(x) = \frac{1}{x^{2} + 8 x + 1}$

$g_8(x) = \frac{1}{\left(x^{2} + 1\right)^{5}}$

$g_9(x) = \frac{1}{x^{8} + x^{2} + 1}$

Exercise 1

y = ((x**3+1)/(x**2+1))**3
Y = ratint(y, x)
display(Latex("Here we consider the function $y=" + sp.latex(y) + "$."))
display(Latex("The integral is $Y=" + sp.latex(Y) + "$."))

Here we consider the function $y=\frac{\left(x^{3} + 1\right)^{3}}{\left(x^{2} + 1\right)^{3}}$.

The integral is $Y=\frac{x^{4}}{4} - \frac{x^{4}}{\left(x^{2} + 1\right)^{2}} + \frac{15 x^{3}}{4 \left(x^{2} + 1\right)^{2}} - \frac{3 x^{2}}{2} - \frac{3 x^{2}}{2 \left(x^{2} + 1\right)^{2}} + 3 x + \frac{13 x}{4 \left(x^{2} + 1\right)^{2}} + 3 \log{\left(\left|{x^{2} + 1}\right| \right)} - \frac{21 \operatorname{atan}{\left(x \right)}}{4}$.

Here we plot all the different terms in $Y$ separately. It is clear that the term $x^4/4$ is larger than all the other terms already when $x=6$, and this effect will become stronger as $x$ increases. To understand this, note that when $x$ is large the $+1$ terms in the definition of $y$ will become negligible, so $y$ will be close to $(x^3/x^2)^3=x^3$, so the main term integral will be $\int x^3\,dx=x^4/4$.

Y_terms = list(Y.args)
Y_terms = [t[1] for t in sorted([(-abs(t.subs({x:5}).evalf()),t) for t in Y_terms])]
xs = np.linspace(0, 6, 100)
fig = go.Figure()
fig.update_layout(width=1000, height=600)
for t in Y_terms:
    t_fun = sp.lambdify(x, t, 'numpy')
    fig.add_trace(go.Scatter(x=xs, y=t_fun(xs), mode='lines', name='$' + sp.latex(t) + '$'))
move_legend(fig, 0.1)
fig.show()

We now plot $Y$ and $x^4/4$ together to see that they are approximately the same.

xs = np.linspace(-9, 9, 100)
Y_fun = sp.lambdify(x, Y, 'numpy')
fig = go.Figure()
fig.update_layout(width=400, height=400)
fig.add_trace(go.Scatter(x=xs, y=Y_fun(xs), mode='lines', name='$Y$'))
fig.add_trace(go.Scatter(x=xs, y=xs**4/4, mode='lines', name='$x^4/4$'))
move_legend(fig, 0.2)
fig.show()

Exercise 2

a, b, c, d = sp.symbols('a b c d', real=True, positive=True)

y = (a*x**2 + b)/(c*x**2 + d)
Y = sp.simplify(sp.integrate(y, x))
display(Latex("Here we consider the function $y=" + sp.latex(y) + "$."))
display(Latex("The integral is $Y=" + sp.latex(Y) + "$."))

Here we consider the function $y=\frac{a x^{2} + b}{c x^{2} + d}$.

The integral is $Y=\frac{a c^{\frac{3}{2}} \sqrt{d} x - c \left(a d - b c\right) \operatorname{atan}{\left(\frac{\sqrt{c} x}{\sqrt{d}} \right)}}{c^{\frac{5}{2}} \sqrt{d}}$.

Y0 = Y.subs({a:7,b:-3,c:2,d:8})
m0 = float((a/c).subs({a:7,b:-3,c:2,d:8}))
display(Y0)
xs = np.linspace(10, 100, 100)
Y0_fun = sp.lambdify(x, Y0, 'numpy')
fig = go.Figure()
fig.update_layout(width=400, height=400)
fig.add_trace(go.Scatter(x=xs, y=Y0_fun(xs), mode='lines', name='$Y_0$'))
fig.add_trace(go.Scatter(x=xs, y=m0*xs, mode='lines', name='$ax/c$'))
move_legend(fig, 0.2)
fig.show()

$\displaystyle \frac{7 x}{2} - \frac{31 \operatorname{atan}{\left(\frac{x}{2} \right)}}{4}$

Exercise 3

u, v, a = sp.symbols('u v a')

y = (8*x+8)/(x**2+2*x+3)
display(y)
display(ratint(y, x))

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

$\displaystyle 4 \log{\left(\left|{x^{2} + 2 x + 3}\right| \right)}$

y = a*(2*x+u)/(x**2+u*x+v)
Y = sp.integrate(y, x)
display(y)
display(Y)

$\displaystyle \frac{a \left(u + 2 x\right)}{u x + v + x^{2}}$

$\displaystyle a \log{\left(u x + v + x^{2} \right)}$

Exercise 3

q = x**2 + 2*x + 3
Delta = q.coeff(x, 1)**2 - 4*q.coeff(x, 2)*q.coeff(x, 0)
print(f"Discriminant = {Delta}")
xs = np.linspace(-2, 1, 100)
ys = sp.lambdify(x, q, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=400, height=400, yaxis_range=[0, 6])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(q) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$\\int ' + sp.latex(1/q) + '\\, dx=' + sp.latex(ratint(1/q, x)) + '$'))

Discriminant = -8

$\int \frac{1}{x^{2} + 2 x + 3}\, dx=\frac{\sqrt{2} \operatorname{atan}{\left(\frac{\sqrt{2} x}{2} + \frac{\sqrt{2}}{2} \right)}}{2}$

q = x**2 + 4*x + 4
Delta = q.coeff(x, 1)**2 - 4*q.coeff(x, 2)*q.coeff(x, 0)
print(f"Discriminant = {Delta}")
xs = np.linspace(-4, 1, 100)
ys = sp.lambdify(x, q, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=400, height=400, yaxis_range=[0, 6])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(q) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$\\int ' + sp.latex(1/q) + '\\, dx=' + sp.latex(ratint(1/q, x)) + '$'))

Discriminant = 0

$\int \frac{1}{x^{2} + 4 x + 4}\, dx=- \frac{1}{x + 2}$

q = x**2 - 5*x + 4
Delta = q.coeff(x, 1)**2 - 4*q.coeff(x, 2)*q.coeff(x, 0)
print(f"Discriminant = {Delta}")
xs = np.linspace(-1, 5, 100)
ys = sp.lambdify(x, q, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=400, height=400, yaxis_range=[-3, 10])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(q) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$\\int ' + sp.latex(1/q) + '\\, dx=' + sp.latex(ratint(1/q, x)) + '$'))

Discriminant = 9

$\int \frac{1}{x^{2} - 5 x + 4}\, dx=\frac{\log{\left(\left|{x - 4}\right| \right)}}{3} - \frac{\log{\left(\left|{x - 1}\right| \right)}}{3}$

Exercise 5

y = g[1]
z = ratint(y, x)
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))
xs = np.linspace(-3, 3, 303)
ix = [(xs < -1),(xs > -1) & (xs < 1),(xs > 1)]
ys = sp.lambdify(x, y, 'numpy')(xs)
zs = sp.lambdify(x, z, 'numpy')(xs)
ys = np.where(np.abs(ys) > 100, np.nan, ys)
zs = np.where(np.abs(zs) > 100, np.nan, zs)
fig = go.Figure()
fig.update_layout(width=800, height=400, yaxis_range=[-10, 10])
for j in range(3):
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=ys[ix[j]], mode='lines', line_color='blue', name='$y$', showlegend=j==0))
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=zs[ix[j]], mode='lines', line_color='red', name='$z$', showlegend=j==0))
fig.add_vline(-1, line_color='green', line_dash='dash')
fig.add_vline( 1, line_color='green', line_dash='dash')

$y=\frac{1}{x^{3} - 1} + \frac{1}{x^{2} - 1} + \frac{1}{x - 1}$

$z=\int y\,dx=\frac{11 \log{\left(\left|{x - 1}\right| \right)}}{6} - \frac{\log{\left(\left|{x + 1}\right| \right)}}{2} - \frac{\log{\left(\left|{x^{2} + x + 1}\right| \right)}}{6} - \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} x}{3} + \frac{\sqrt{3}}{3} \right)}}{3}$

y = g[2]
z = ratint(y, x)
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))
xs = np.linspace(-2, 5, 1000)
ix = [(xs < 1),(xs > 1) & (xs < 2),(xs > 2) & (xs < 3),(xs > 3)]
ys = sp.lambdify(x, y, 'numpy')(xs)
zs = sp.lambdify(x, z, 'numpy')(xs)
ys = np.where(np.abs(ys) > 100, np.nan, ys)
zs = np.where(np.abs(zs) > 100, np.nan, zs)
fig = go.Figure()
fig.update_layout(width=800, height=400, yaxis_range=[-50, 50])
for j in range(4):
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=ys[ix[j]], mode='lines', line_color='blue', name='$y$', showlegend=j==0))
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=zs[ix[j]], mode='lines', line_color='red', name='$z$', showlegend=j==0))
for x0 in [1, 2, 3]:
    fig.add_vline(x0, line_color='green', line_dash='dash')

fig.show()

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

$z=\int y\,dx=- \frac{1}{x - 1} - \frac{1}{2 \left(x - 2\right)^{2}} - \frac{1}{3 \left(x - 3\right)^{3}}$

y = g[3]
z = ratint(y, x)
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))
xs = np.linspace(-5, 5, 1000)
ix = [(xs < -1),(xs > -1)]
ys = sp.lambdify(x, y, 'numpy')(xs)
zs = sp.lambdify(x, z, 'numpy')(xs)
ys = np.where(np.abs(ys) > 100, np.nan, ys)
zs = np.where(np.abs(zs) > 100, np.nan, zs)
fig = go.Figure()
fig.update_layout(width=800, height=400, yaxis_range=[-10, 10])
for j in range(2):
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=ys[ix[j]], mode='lines', line_color='blue', name='$y$', showlegend=j==0))
    fig.add_trace(go.Scatter(x=xs[ix[j]], y=zs[ix[j]], mode='lines', line_color='red', name='$z$', showlegend=j==0))
for x0 in [-1]:
    fig.add_vline(x0, line_color='green', line_dash='dash')

fig.show()

$y=\frac{1}{x^{3} + 1} + \frac{1}{x^{2} + 1} + \frac{1}{x + 1}$

$z=\int y\,dx=\frac{4 \log{\left(\left|{x + 1}\right| \right)}}{3} - \frac{\log{\left(\left|{x^{2} - x + 1}\right| \right)}}{6} + \operatorname{atan}{\left(x \right)} + \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} x}{3} - \frac{\sqrt{3}}{3} \right)}}{3}$

Exercise 6

y = g[2]
d = sp.denom(sp.factor(y))
display(Latex(f"For $y={sp.latex(y)}$ the denominator is $d={sp.latex(d)}$"))
xs = np.linspace(0, 4, 500)
ys = sp.lambdify(x, d, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=600, height=400, yaxis_range=[-0.3, 0.3])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(d) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))

For $y=\frac{1}{\left(x - 1\right)^{2}} + \frac{1}{\left(x - 2\right)^{3}} + \frac{1}{\left(x - 3\right)^{4}}$ the denominator is $d=\left(x - 3\right)^{4} \left(x - 2\right)^{3} \left(x - 1\right)^{2}$

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

$z=\int y\,dx=\frac{4 \log{\left(\left|{x + 1}\right| \right)}}{3} - \frac{\log{\left(\left|{x^{2} - x + 1}\right| \right)}}{6} + \operatorname{atan}{\left(x \right)} + \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} x}{3} - \frac{\sqrt{3}}{3} \right)}}{3}$

y = g[3]
d = sp.denom(sp.factor(y))
display(Latex(f"For $y={sp.latex(y)}$ the denominator is $d={sp.latex(d)}$"))
xs = np.linspace(-2, 2, 500)
ys = sp.lambdify(x, d, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=600, height=400, yaxis_range=[-10, 10])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(d) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))

For $y=\frac{1}{x^{3} + 1} + \frac{1}{x^{2} + 1} + \frac{1}{x + 1}$ the denominator is $d=\left(x + 1\right) \left(x^{2} + 1\right) \left(x^{2} - x + 1\right)$

$y=\frac{1}{x^{3} + 1} + \frac{1}{x^{2} + 1} + \frac{1}{x + 1}$

$z=\int y\,dx=\frac{4 \log{\left(\left|{x + 1}\right| \right)}}{3} - \frac{\log{\left(\left|{x^{2} - x + 1}\right| \right)}}{6} + \operatorname{atan}{\left(x \right)} + \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} x}{3} - \frac{\sqrt{3}}{3} \right)}}{3}$

y = g[5]
d = sp.denom(sp.factor(y))
display(Latex(f"For $y={sp.latex(y)}$ the denominator is $d={sp.latex(d)}$"))
xs = np.linspace(-4, 1, 500)
ys = sp.lambdify(x, d, 'numpy')(xs)
fig = go.Figure()
fig.update_layout(width=600, height=400, yaxis_range=[-2, 15])
fig.add_trace(go.Scatter(x=xs, y=ys, mode='lines', name='$' + sp.latex(d) + '$', showlegend=True))
move_legend(fig, 0.2)
fig.show()
display(Latex('$y=' + sp.latex(y) + '$'))
display(Latex('$z=\\int y\\,dx=' + sp.latex(z) + '$'))

For $y=\frac{x \left(x + 2\right) \left(x + 4\right)}{\left(x + 1\right)^{2} \left(x + 3\right)^{2}}$ the denominator is $d=\left(x + 1\right)^{2} \left(x + 3\right)^{2}$

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

$z=\int y\,dx=\frac{4 \log{\left(\left|{x + 1}\right| \right)}}{3} - \frac{\log{\left(\left|{x^{2} - x + 1}\right| \right)}}{6} + \operatorname{atan}{\left(x \right)} + \frac{\sqrt{3} \operatorname{atan}{\left(\frac{2 \sqrt{3} x}{3} - \frac{\sqrt{3}}{3} \right)}}{3}$

def randpoly(x, degree=3, coef_range=(-10, 10)):
    N = degree + 1
    a = np.random.randint(*coef_range, N)
    f = sp.Add(*[a[i] * x **i for i in range(N)])
    return f    

x = sp.symbols('x', real=True)

xs = np.pi * np.linspace(-5, 5, 10) + np.exp(1)
f = randpoly(x, 3) / randpoly(x, 4)
F = ratint(f, x)
err0 = sp.simplify(sp.diff(F,x) - f)
err1 = [np.abs(err0.subs(x, x0).evalf()) for x0 in xs]
err2 = np.max(err1)
display(f)
display(F)
print(err2)

$\displaystyle \frac{- 5 x^{3} + x^{2} + x}{- 9 x^{4} + 5 x^{3} + 6 x^{2} + 9 x + 9}$

$\displaystyle 0.201469843915853 \log{\left(\left|{x - 1.60317241722703}\right| \right)} + 0.0848043889509342 \log{\left(\left|{x + 0.777136137892559}\right| \right)} + 0.134640661344384 \log{\left(\left|{x^{2} + 0.270480723778914 x + 0.802643444735052}\right| \right)} - 0.0261255625968842 \operatorname{atan}{\left(1.1291303566521 x + 0.152703998054001 \right)}$

3.12286456336813e-17