MAS2008 Scientific Computing: Lab 5
Integration, root finding


The following notebooks and videos are relevant for this lab:
Borwein integrals View Download
Numerical solution of the Laplace equation View Download
There are also two interactive demonstrations: one for quadrature methods and one for Laplace's equation.

Task 1: Integration over the unit disc

This task is about numerical evaluation of $\iint_D f(x,y) dxdy$, where $f(x,y)$ is a function defined on the unit disc $D$ where $x^2+y^2\le 1$.

First, define $J(i,j)=0$ if $i$ or $j$ is odd, and \[ J(2k,2l) = J(2l,2k) = \frac{\pi}{4^{k+l}(k+l+1)}\sum_{m=l-k}^{l+k} (-1)^{l+m} \binom{2l}{m} \binom{2k}{k+l-m} \] if $k\leq l$. It can be shown that $J(i,j)$ is the exact value of the integral of $x^i y^j$ over $D$. Define a function disc_int_monomial(i, j) that returns $J(i,j)$. You should just use the formula above, and not try to compute the integral directly. You might need to search the web or ask an AI assistant how to deal with binomial coefficients in Python.

For numerical integration, we will use various rules of the form \[ \iint_D f(x,y)dx\,dy \approx \sum_{i=1}^m w_i f(x_i,y_i) \] where $(x_i,y_i)$ are points in $D$ and $w_i$ are weights. Usually we are interested in symmetric rules, so if we have a point $(x,y)$ with weight $w$ then we also have $(-x,y)$, $(x,-y)$, $(-x,-y)$, $(y,x)$, $(-y,x)$, $(y,-x)$, and $(-y,-x)$ with the same weight. For a typical point $(x,y)$, this gives 8 points in total. However, in special cases some of these points may coincide: Define a function extend_disc_rule(u0) as follows:

Now define a function show_disc_rule(u) that produces a plot of a rule u consisting of rows [x,y,w,t] as produced by extend_disc_rule().

Now define a function disc_int(f,u) that takes a function f(x,y) and a rule u as produced by extend_disc_rule(), and returns the corresponding numerical approximation to $\iint_D f(x,y) dxdy$. You should use the formula \[ \iint_D f(x,y) dxdy \approx \sum_{i=1}^m w_i f(x_i,y_i) \] where $(x_i,y_i)$ are the points in the rule and $w_i$ are the weights.

Now try disc_int(lambda x,y: x**2 * y**4, kim_song_rule) and compare the answer with disc_int_monomial(2,4), then repeat for various other exponents instead of 2 and 4. You should find that the answers are extremely close when the sum of the exponents is less than 10, but somewhat worse when the exponents are higher.

Task 2: Roots of the Wilkinson polynomial

Here we consider Wilkinson's polynomial \[ p_n(x) = (x-1)(x-2)\dotsb(x-n) = \prod_{i=1}^n(x-i), \] together with the perturbed polynomial $p_{n,\epsilon}(x)=p_n(x)+\epsilon\,x^{n-1}$. Clearly, the roots of $p_n(x)$ are $x=1,\dotsc,n$, so if $\epsilon$ is sufficiently small, then the roots of $p_{n,\epsilon}$ will be close to $x=1,\dotsc,n$. However, this is a good example of a problem where numerical algorithms struggle to locate the roots accurately.

Polynomials can be represented in Python as numpy arrays of coefficients, for example $2x^3-3x^2+4x-5$ corresponds to np.array([2,-3,4,-5]) and $x^4-1$ corresponds to np.array([1,0,0,0,-1]). If $f(x)$ corresponds to $u$ and $g(x)$ corresponds to $v$ then the product $f(x)g(x)$ corresponds to np.convolve(u,v). The function np.roots(u) returns the roots of the polynomial corresponding to $u$.

Using this, write a function wilkinson_poly(n,epsilon) that returns the coefficients of the polynomial $p_{n,\epsilon}(x)$. In this function and all the functions discussed later, the parameter $\epsilon$ should default to $0$ if not specified. Then write a function wilkinson_roots(n,epsilon) that returns the roots of $p_{n,\epsilon}(x)$. Except when $\epsilon$ is tiny, some of the roots will be complex. Write a function wilkinson_roots_argand_plot(n,epsilon) that plots the roots of $p_{n,\epsilon}(x)$ in the complex plane.

Now fix $n$ and consider how the roots change as $\epsilon$ varies. Define a function wilkinson_roots_real_part_plot(n) that plots the points $(t,x)$ where $t$ varies from $-10$ to $1$ in 100 steps and $x$ runs over the real parts of the roots of $p_{n,\exp(t)}(x)$. Also do the corresponding thing for the imaginary parts of the roots.