I am going to discuss my reasons for wanting this first, as this may in fact not be what I am looking for.

My reason for asking this that I have finished writing a piece of code that solves, $-\nabla \cdot(a \nabla u) = f$ where $a$ is discontinuous, this is the result of the code

enter image description here

we can see oscillations about the line of discontinuity, so I have been told that I have to enrich the approximation. After some manipulation, the meat of the problem is dealing with this integral:

$$\int_{\Omega}a\nabla u \cdot \nabla v dx dy$$

On enriching my approximation.

$$u(x,y) = \sum_{i = 0}^N\sum_{j = 0}^Nu_{ij}h_i(x)h_j(y) + \sum_{k = 1}^N\sum_{l = 1}^N\alpha_{kl}\tilde{h_k(x)}\tilde{h_l(y)}\Phi(x,y)$$

On substituting this into the integral

$$\int_{\Omega}a\nabla u \cdot \nabla v dx dy $$ $$\approx w_s\sum_{n = 0}^N \sum_{i = 0}^N w_n h_i'(x_n)h_r'(x_n)u_{is}, \forall s\in[0,N]$$

$$+ w_r\sum_{m = 0}^N \sum_{j = 0}^N w_m h_j'(y_m)h_s'(y_m)u_{rj}, \forall r\in[0,N]$$

$$+ \color{red}{ w_s\sum_{n = 1}^{2N + 2} a_1(x_n,y_s)\sum_{k = 1}^{N-1} w_n \tilde{h_k}'(x_n)h_r'(x_n)\alpha_{ks}\Phi_{1,x}(x_n,y_s), \forall s\in[1,N]}$$

$$+ w_r\sum_{m = 0}^{2N + 2} a_1(x_r,y_m) \sum_{l = 0}^{N-1} w_m \tilde{h_l}'(y_m)h_s'(y_m)\alpha_{rl}\Phi_{1,y}(x_r,y_m), \forall r\in[1,N]$$

$$+ \color{red}{w_s\sum_{n = 1}^{2N + 2} a_2(x_n,y_s)\sum_{k = 1}^{N-1} w_n \tilde{h_k}'(x_n)h_r'(x_n)\alpha_{ks}\Phi_{2,x}(x_n,y_s), \forall s\in[1,N]}$$

$$+ w_r\sum_{m = 0}^{2N + 2} a_2(x_r,y_m) \sum_{l = 0}^{N-1} w_m \tilde{h_l}'(y_m)h_s'(y_m)\alpha_{rl}\Phi_{2,y}(x_r,y_m), \forall r\in[1,N] \textbf{(1)}$$

How do I form the Chebyshev differentiation matrix in MATLAB?

The first answer in this link deals with the definition of the $h_i'(x_j)$ functions, or the construction of the square Legendre pseudospectral differentiation matrix. The $\tilde{h_i}'(x_j)$ have a slightly different definition.

$$h_i'(x) = \frac{(1-x^2)(L_N'(x))}{N(N+1)L_N'(x_i)(x - x_i)}$$ $$\tilde{h_i}'(x) = \frac{(1-x_i^2)(L_N'(x))}{N(N+1)L_N'(x_i)(x - x_i)}$$

The construction of the corresponding square Legendre pseudospectral differentiation matrix for each of the above is given in code as


if N<3

%N1xN1 zero matrix
%eps= epsilon!
while max(abs(x-xold))>eps

    %set first collumn entries to 1 (P(:,1) = 1);
    %set second collumn entries to x (P(:,2) = x);
    % and so on, the collumns of this matrix are the legendre polynomials
    for k=2:N
        P(:,k+1)=( (2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1) )/k;
        %%Bonnets formula ;)
    % Roots of (1-x^2)L'_N
    x=xold-( x.*P(:,N1)-P(:,N) )./( N1*P(:,N1));

%---Legendre pseduspectral differentiation matrix------------page 124-----
% of the book 

x = -x;
for k=1:N+1;
    for j=1:N+1;
D(1,1)=-N*(N+1)/4 ;           

C1 = -(1-x(1:N1).^2);
C2 = N*N1.*P(:,N1);
C = C1+C2;
% L1 = P' L2 = P''...
% LN(i,N) is the Nth derivative of the Nth legendre polynomial at x(i)

% The last collumn in each is only really relevant.
%The coll's before will give a wierd result, this is easily fixed 
%but we dont use them so....

L1 = zeros(N1,N1);
L2 = zeros(N1,N1);
L3 = zeros(N1,N1);
L1(:,1) = 0;
L1(:,2) = 1;

for m = 2:N;
    for k = 1:N1;
        L1(k,m+1) = ((2*m + 1)*(P(k,m) + x(k)*L1(k,m)) - m*L1(k,m-1))/(m+1); 
L1(1,N1) = N*N1/2;
L1(N1,N1) = (N*N1/2)*(-1)^(N-1);

L2(:,1) = 0;
L2(:,2) = 0; 
for m = 2:N;
    for k = 1:N1;
        L2(k,m+1) = ((2*m + 1).*(2*L1(k,m) + x(k).*L2(k,m)) - m.*L2(k,m-1))/(m+1); 
L2(1,N1) = (N*N1 - 2)*(N*N1)^2/8;
L2(N1,N1) = (N*N1 - 2)*((N*N1)^2/8)*(-1)^N;

L3(:,1) = 0;
L3(:,2) = 0;
for m = 2:N;
    for k = 1:N1;
        L3(k,m+1) = ((2*m + 1).*(3*L2(k,m) + x(k).*L3(k,m)) - m.*L3(k,m-1))/(m+1); 

L3(1,N1) = (N*N1 - 6)*L2(1,N1)/6;
L3(N1,N1) = ((N*N1 - 6)*L2(1,N1)/6)*(-1)^N;
%this is the matrix that corresponds to the derivatives of h tilde.
DT = zeros(N1,N1);
for l = 2:N;
    for k = 2:N;
        if(k == l)
            DT(l,k) = C(k)*(L3(k,N1)/2);
            DT(l,k) = C(k)*(L2(l,N1)/(x(l) - x(k)));
for k = 2:N;
    DT(1,k) = C(1)*((1 - x(k))*(L2(1,N1) - L1(1,N1))/(1-x(k))^2); 
    DT(N1,k) = C(N1)*((1 - x(k))*(L2(N1,N1) - L1(N1,N1))/(1-x(k))^2); 

The reason for the 4 extra terms is that I have split the region into two, either side of the discontinuity. According to my professor I have to use a technique called over-integration. I am going to need to define another D matrix such that is it a $(2N+2) \times (N+1)$ matrix and similarly for DT(ilde).

So my question is

How do I form the non-square Legendre pseudospectral differentiation matrix?, if such a thing exists. If not, how do I deal with the over integration parts of the expression?

  • 1
    $\begingroup$ Apologies for not exactly answering your question, but have you considered using e.g. discontinuous Galerkin methods? They are specifically developed to deal with coefficient jumps and are quite simple to implement. If your objective is not the use of a specific method it may be simpler than enriching the standard space. $\endgroup$
    – moyner
    Sep 13, 2016 at 8:13
  • $\begingroup$ You might be interested to check your solution against a straight forward implementation within the chebfun framework. In regard of spectral methods I can also recommend the literature referenced there. $\endgroup$
    – Bort
    Sep 14, 2016 at 12:10
  • $\begingroup$ I got so stuck on this, I had given up. I think the main issue was that I cant seem to define the matrix described below in MATLAB. I will look through this though, thanks for the suggestion. I think the implementation of the D tilde matrix is also incorrect. Just a mess really. $\endgroup$
    – user21030
    Sep 14, 2016 at 14:33

1 Answer 1


For the purpose of solving a differential equation, the differentiation matrix must be square (and invertible). You have added additional degrees of freedom -- thus increasing the number of columns -- so you also need to add the same number of (independent) conditions by increasing the number of rows.

If you have selected your enriching functions properly, this can be done by using these functions as additional test functions, i.e., inserting $$v_{N+k,N+l}(x,y) = \tilde h_k(x)\tilde h_l(y)\Phi(x,y)$$ and integrating. Since LGL-pseudospectral methods use nodal basis functions based on quadrature points, this can be seen as adding more quadrature points for each element; in particular, since the enriching functions do not raise the maximal degree, more than you need for exact integration of each basis function. Hence, this is sometimes referred to as over-integration; see, e.g., this reference.

  • $\begingroup$ Just so I am clear before I attempt to code this. I have edited the question above to highlight two terms of expression $\textbf(1)$, so these two terms are $\int_{\Omega_2} \frac{\partial \sum_{k = 1}^N\sum_{l = 1}^N\alpha_{kl}\tilde{h_k(x)}\tilde{h_l(y)}\Phi(x,y)}{\partial x} \frac{\partial h_r(x) h_s(y)}{\partial x}$. they are split into two. So regarding the differentiation matrices we define two square differentiation matrices $DT_M$ $D_M$ each of size $(2N+2) \times (2N+2) $ and when programming the first red term $k,r \in [0,N+1]$ then in the second $k,r \in [N+1,2N+2]$? $\endgroup$
    – user21030
    Sep 13, 2016 at 11:06
  • $\begingroup$ No, you always need all terms when substituting $u$. Only for the test functions $v$ it suffices to fix all but one of the $u_{ij},\alpha_{kl}$ at $0$ (and the remaining one at $1$). I'd suggest not thinking of this as splitting the matrix, but just pretend you now have $2(N+1)$ basis functions instead of $N+1$. $\endgroup$ Sep 13, 2016 at 11:08
  • $\begingroup$ I am sorry about this, I am trying to understand this. So I am correct in saying that I need to define a $(2N+2) \times (2N+2)$ version of the differentiation matrix? ($D_{M = 2N+2}$)This I must have to evaluate $h'$ at $n\in[1,2N+2]$ points( the rows). from the ranges I have defined above I would only be using upto the $(N+1)th$ column of $D_{M}$. Would this be correct? $\endgroup$
    – user21030
    Sep 13, 2016 at 12:40
  • $\begingroup$ Yes, if you mean what I think you mean. In the end, you need to end up with a $(2N+2)\times (2N+2)$ matrix -- of which the upper left $(N+1)\times (N+1)$ block is your original, unenriched differentiation matrix. Evaluating the new basis functions at the old points gives you another $(N+1)\times (N+1)$ upper right block. Now you need a lower left (old basis function at new evaluation points) and a lower right (new basis functions at new points) block of the same size to complete the square. $\endgroup$ Sep 13, 2016 at 12:54
  • $\begingroup$ I am confused what you mean by old points and new points. Really confused. So the matrix you described above is a matrix that contains the entries of $h'(x_n)$ and $\tilde{h}'(x_n)$? or are you describing the stiffness matrix? It has been a few days, I am sorry. $\endgroup$
    – user21030
    Sep 16, 2016 at 15:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.