Solving double integral numerically in matlab

Question

In the paper "The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator". Where the author has solved a fractional laplacian equation on bounded domain as a non-local diffusion equation.

I am trying to implement the finite element approximation of the one dimensional problem(please refer to page 14 of the above mentioned paper) in matlab.

I am using the following definition of $\phi_k$ as it is mentioned in the paper that $\phi$ is a $hat\;function$ \begin{equation} \phi_{k}(x)=\begin{cases} {x-x_{k-1} \over x_k\,-x_{k-1}} & \mbox{ if } x \in [x_{k-1},x_k], \\ {x_{k+1}\,-x \over x_{k+1}\,-x_k} & \mbox{ if } x \in [x_k,x_{k+1}], \\ 0 & \mbox{ otherwise},\end{cases} \end{equation}

$\Omega=(-1,1)$ and $\Omega_I=(-1-\lambda,-1) \cup (1,1+\lambda)$ so that $\Omega\cup\Omega_I=(-1-\lambda,1+\lambda)$

For the integers K,N we define the partition of $\overline{\Omega\cup\Omega_I}=[-1-\lambda,1+\lambda]$ as,

\begin{equation} -1-\lambda=x_{-K}<...

Finally the equations that we have to solve to get the solution $\tilde{u_N}=\sum_{i=-K}^{K+N}U_j\phi_j(x)$ for some coefficients $U_j$ is:

enter image description here

Where $i=1,...,N-1$.

I need pointers in order to simplify and solve the LHS double integral in matlab.It is written in the paper(page 15) that I should use four point gauss quadrature for inner integral and quadgk.m function for outer integral, but since the limits of the inner integral are in terms of x how can I apply four point gauss quadrature on it??.Any help will be appreciated. Thanks.

You can find the original question here.(Since SO does not support Latex)

Suggest math.stackexchange.com instead (yah, I know, it's programming, but you may have more luck there.) — Jason S
I suggest you split this question up in two parts: 1) ask for the simplification on mathoverflow, 2) ask for an implementation here on SO. On SO, leave out all the math except the double integral (I put a picture in for you). On mathoverflow, leave out the request for an implementation in MATLAB. — Rody Oldenhuis

Rody Oldenhuis Rody Oldenhuis · Accepted Answer · 2013-10-21T07:17:01

For a first stab at the problem, take a look at dblquad and/or quad2d.

In the end, you'll want custom quadrature methods, so you should do something like the following:

% The integrand is of course a function of both x and y
integrand = @(x,y) (phi_j(y) - phi_j(x))*(phi_i(y) - phi_i(x))/abs(y-x)^(2*s+1);

% The inner integral is a function of x, and integrates over y
inner = @(x) quadgk(@(y)integrand(x,y), x-lambda, x+lambda);

% The inner integral is integrated over x to yield the value of the double integral 
dblIntegral = quadgk(inner, -(1+lambda), 1+lambda)

where I've used quadgk twice, but you can replace by any other (custom) quadrature method you please.

By the way -- what is the reason for the authors to suggest a (non-adaptive) 4-point Gauss method? That way, you have no estimation of (and/or control over) the errors made in the inner integral...

Solving double integral numerically in matlab

2 Answers