An important operation in numerical analysis is to calculate the integrals of functions. In general, $$ \int_{-1}^{1} u(x)\:dx \approx \sum_{i=1}^N w_i u(x_i) $$ In this expression, $w_i$ is a weight for node $x_i$. For different numerical integration procedures, different weights and nodes are used.
If polynomials of order $m$ are to be integrated exactly, initially, since $m+1$ coefficients are necessary to define the polynomial, there should be $m+1$ arbitrary nodes and weights. If, on the other hand, the nodes can be fixed beforehand, it is possible to do better.
Often a weighed integral is useful:
$$ \int_{-1}^{1} (1-x)^\alpha(1+x)^\beta u(x)\:dx = \sum_{i=1}^Q w_i^{\alpha,\beta} u(x^{\alpha,\beta}_i) + R(u) $$where R(u) is a remainder which can be zero if $u(x) \in \mathcal{P}_{2Q-k}([-1,1])$
Depending on the number of nodes that are free to float (assume any value necessary to the procedure) different integration schemes are possible:
gj.grjm when -1 is included or grjp when 1 is included.glj.Nodes emplyed in the quadrature can be calculated using the following functions:
zgj(Q, a, b) for Gauss-Jacobi with $Q$ nós, weights $\alpha=$a and $\beta=$b.zgj(Q, a, b) for Gauss-Lobatto-Jacobi.zgrjm(Q, a, b) for Gauss-Radau-Jacobi including the end -1.zgrjp(Q, a, b) for Gauss-Radau-Jacobi including the end 1.
In [1]:
using Jacobi
Q = 5
alpha=0.3
beta = 0.8
zgj(Q, alpha, beta) # Gauss-Jacobi
Out[1]:
In [2]:
zgrjm(Q, alpha, beta) # Gauss-Radau-Jacobi, -1
Out[2]:
In [3]:
zgrjp(Q, alpha, beta) # Gauss-Radau-Jacobi, 1
Out[3]:
In [4]:
zglj(Q, alpha, beta) # Gauss-Lobatto-Jacobi
Out[4]:
The default case is to calculate the zeros using Float64, if some other type is necessary, just add the type as a fourth argument. The only restriction is that this type should be a FloatingPoint since iteration is used.
In [5]:
zglj(5, 0, 0, Float32)
Out[5]:
In [6]:
zgj(5, 0, 0, BigFloat)
Out[6]:
In [7]:
a = 0
b = 0
z = zgj(5, a, b)
w = wgj(z, a, b)
Out[7]:
To compute an integral,
In [8]:
using Polynomials
u = Poly([10,9,8,7,6,5,4,3,2,1])
Out[8]:
In [9]:
Iu = polyint(u)
Out[9]:
In [10]:
I = polyval(Iu, 1.0) - polyval(Iu, -1.0)
Out[10]:
In [11]:
ui = polyval(u, z)
Out[11]:
In [12]:
Igauss = sum(w.*ui)
Out[12]:
In [13]:
err = Igauss - I
Out[13]:
Since the quadrature had Q=5 nodes and a Gauss-Jacobi scheme was used, the exact value of the integral was possible. As another example let's use the weights $\alpha=1$ and $\beta=1$. For integer weights, the integrand is stilla polynomial:
In [14]:
p = Poly([1,-1]) * Poly([1,1]) * u
Out[14]:
Using the same weights,
In [15]:
Ip = polyint(p)
I2 = polyval(Ip, 1.0) - polyval(Ip, -1.0)
Igauss2 = sum(w.*polyval(p, z))
I2, I2-Igauss
Out[15]:
With new weights:
In [16]:
Q = 5
z2 = zgj(Q, 1, 1)
w2 = wgj(z2, 1, 1)
ui = polyval(u, z2)
Igab = sum(w2.*ui)
I2 - Igab
Out[16]:
Which is zero, except for some floting point error!
The nodes of the quadrature can be used for polynomial interpolation using Lagrangean interpolants. If a function us evaluated at the nodes of the quadrature,
$$ u_i = u(x_i^{\alpha,\beta}) $$The corresponding Lagrangean interpolant is $h_i(x)$:
$$ h_i(x) = \prod_{k=1,k\ne i}^Q \frac{x - x_k^{\alpha,\beta}}{x_i^{\alpha,\beta} - x_k^{\alpha,\beta}} $$This can be computed using function lagrange(i, x, z) where z is an array of interpolation nodes
In [17]:
using PyPlot
x = linspace(-1,1, 401)
np = length(x)
Q = 5
z = zglj(Q)
h = zeros(np, Q)
for i = 1:Q
for k = 1:np
h[k,i] = lagrange(i, x[k], z)
end
plot(x, h[:,i])
axvline(x=z[i], color="black")
end
axhline(y=0, color="black");
This was calculated using Gauss-Lobatto-Jacobi nodes. If equally spaced nodes were used things are not as pretty:
In [18]:
z = linspace(-1, 1, Q)
h = zeros(np, Q)
for i = 1:Q
for k = 1:np
h[k,i] = lagrange(i, x[k], z)
end
plot(x, h[:,i])
axvline(x=z[i], color="black")
end
axhline(y=0, color="black");
It still works well but as the number of nodes increases, it will get worse.
In [19]:
Q = 9
x = linspace(-1,1,401)
z = zglj(Q)
u = sin(2pi*z)
ue = sin(2pi*x)
imat = interp_mat(x, z)
ui = imat * u
plot(x, ue, "r", z, u, "rs", x, ui, "b--");
In [20]:
x = linspace(-1,1, 401)
Q = 10
z = zglj(Q)
u = z .* cos(2pi*z)
Dmat = dglj(z)
du = Dmat * u
du_exact = cos(2pi*x) - 2pi*x .* sin(2pi*x)
plot(x, du_exact, "r", z, du, "rs");
In [21]:
q = Quadrature(GLJ, 4, 1, 1)
Out[21]:
In [22]:
qzeros(q)
Out[22]:
In [23]:
qweights(q)
Out[23]:
In [24]:
qdiff(q)
Out[24]:
In [25]:
qalpha(q), qbeta(q)
Out[25]:
In [ ]:
The differential equation is
$$ \frac{d^2u}{dx^2} + u = \sin 2\pi x $$To make things simpler, we will use the following Neumann BCs:
$$ \left.\frac{du}{dx}\right|_{x=-1} = \left.\frac{du}{dx}\right|_{x=1} = 0 $$The weak form, for appropriate function spaces is:
$$ \int_{-1}^1 \frac{dv}{dx}\frac{du}{dx}\:dx - \left.v\frac{du}{dx}\right|_{-1}^1 - \int_{-1}^1 uv\:dx = -\int_{-1}^1 vf\:dx $$Using Lagrange interpolation at the nodes of the Gauss-Lobatto-Jacobi quadrature, we can solve this equation:
$$ u(x) = \sum_{k=1}^Q u_kh_k(x) \qquad v(x) = \sum_{i=1}^Q u_ih_i(x) $$Substituting in the weak form,
$$ \left([L] - [M]\right)\{u\} = -[M]\{f\} $$where
$$ M_{ik} = \int_{-1}^1 h_kh_i\:dx $$Usando a propria quadratura utilizada para definir os interpoladores $h_i(x)$ para se integrar esta expressão se obtém uma aproximação da integral mas próxima o suficiente se Q for grande. Com isso, esta matriz é diagonal com cada elemento da diagoal igual ao peso da quadratura: $$ M_{ik} = w_i\delta_{ik} $$ Ja a outra matrix é um pouco mais complicado: $$ L_{ik} = \int_{-1}^1 \frac{dh_i}{dx}\frac{dh_k}{dx}\:dx = \sum_{\alpha=1}^Q w_\alpha D_{\alpha i} D_{\alpha k} $$
In [26]:
function solve(Q, fun, x)
z = zglj(Q)
w = wglj(z)
D = dglj(z)
imat = interp_mat(x, z)
f = fun(z)
L = zeros(Q,Q)
for i=1:Q, k=1:Q
ll = 0.0
for a = 1:Q
ll = ll + w[a] * D[a,i] * D[a, k]
end
L[k,i] = ll
end
fr = zeros(Q)
for i = 1:Q
L[i,i] = L[i,i] - w[i]
fr[i] = - f[i]*w[i]
end
u = L\fr
ue = imat * u
return ue
end
Out[26]:
In [27]:
Q = 3:20
nq = length(Q)
x = linspace(-1, 1, 101)
np = length(x)
u = zeros(np, nq)
for i = 1:nq
u[:,i] = solve(Q[i], x->sin(2pi*x), x)
end
for i = 1:(nq-1)
plot(x, u[:,i])
end
plot(x, u[:,nq], "--", linewidth=4, color="black");
In [28]:
# Cálculo do erro:
ue = u[:,nq]
uerr = zeros(nq-1)
for i = 1:(nq-1)
uerr[i] = maxabs(u[:,i] - ue)
end
plot(log10(Q[1:(nq-1)]), log10(uerr), "rs")
xlabel("log(Q)")
ylabel("log(error)");
In [ ]: