MAS2008 Scientific Computing: Lab 5 Integration, root finding

Instructions

The following notebooks and videos are relevant for this lab:
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:
• If $x=y=0$ then all 8 points in the list are equal to $(0,0)$. We say that this point has type 0.
• If $y=0$ but $x\neq 0$ then we only have four points $(\pm x,0)$ and $(0,\pm x)$. We say that points like this have type 1.
• If $x=y\neq 0$ then we only have four points $(\pm x,\pm x)$. We say that points like this have type 2.
• If $|x|\neq |y|$ and $x,y\neq 0$ then we have 8 distinct points. We say that points like this have type 3.
Define a function extend_disc_rule(u0) as follows:
• You can assume that u0 will be a numpy array of shape $(n,3)$ for some $n\geq 1$. Each row will have the form [x,y,w], where $(x,y)$ has type 0, 1, 2 or 3 as described above, and $w$ is a weight.
• Your function should return a numpy array of shape $(m,4)$ for some $m\geq n$. Each row will have the form [x,y,w,t], where $w$ is a weight and $t$ is the type of $(x,y)$. Each row of type 3 in the input should give rise to 8 rows of type 3 in the output, each row of type 1 or 2 in the input should give rise to 4 rows in the output, and each row of type 0 in the input should give rise to a single row in the output.
• You should work out the case [[0,0,0.2],[0.5,0,0.05],[0.625,0.625,0.05],[0.875,0.125,0.05]] by hand, and check that your function gives the expected output.
• You can also cut and paste the definition of kim_song_rule0 from this file, and then put kim_song_rule=extend_disc_rule(kim_song_rule0). You should find that kim_song_rule has 57 rows. (This comes from the paper A survey of known and new cubature formulas for the unit disk by Cools and Kim, which cites an earlier paper by Kim and Song.)

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().
• The plot should be square with the same scale for the two axes, but no axes should be visible.
• The plot should show the unit circle. You can either use formula $(\cos(t),\sin(t))$ for $t\in[0,2\pi]$, or read the documentation for matplotlib.patches.Circle.
• For each point $(x,y,w,t)$ in the rule, you should draw a filled circle centred at $(x,y)$. The radius should be proportional to $w$. If you use plt.scatter then the size must be specified in typographic points, which are small: a size of $400w$ points is about right. Points of type 0,1,2 or 3 should be black, red, green and blue respectively.

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.