Borwein integrals | View | Download | |
Numerical solution of the Laplace equation | View | Download |
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.
extend_disc_rule(u0)
as follows:
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.
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.
[[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.
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.)
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()
.
matplotlib.patches.Circle
.
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.
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.
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.
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$.
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.
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.