Integration, root finding

Borwein integrals | View | Download | |

Numerical solution of the Laplace equation | View | Download |

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.

`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.
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.