National Taiwan University
Check your differentiation matrix!
Note the bottom left part!
Also check your mesh don't cover both end points!
You can choose one of the ways below.
Given a sequence $x_k$ which converges to $L$ and denote the error $(L - x_k)$ by $e_k$.
Then the rate of convergence is defined by
$$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|}=c.$$
And can be categorized into 3 cases
Linearly $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|}=c \in (0,1)$$
Superlinearly $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|}=0$$
Sublinearly $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|}=1$$
$$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|^n} < M \in (0, \infty)$$
And a sequence is called Converge Logarithmically if
it converge sublinearly $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|} = 1$$
and $$\lim_{k\to\infty}\frac{|e_{k+2}-e_{k+1}|}{|e_{k+1}-e_k|} = 1$$
Also you can distinguish them roughly by the plot
From left to right: linearly, linearly, superlinearly, sublinearly (Ref. Wikipedia)
You can evaluate the rate of convergence by the following sample code.
(Here we use 2-norm as example)
So the rate of convergence is about 0.35356615, which is about $2^{-3/2}$
So... the right answer is "it converges linearly."
(NOT superlinearly, sublinearly, logarithmically...)
We know for the one-sided 3-point differentiation approximation
$$f'(x) \approx g(x) \triangleq \alpha_1 f(x) + \alpha_2 f(x+h) + \alpha_3 f(x+2h)$$
can be derived by using the Taylor expansion and then solving the linear system $$\begin{aligned} f(x) &= f(x)\\ f(x+h) &= f(x) + hf'(x) + \frac{h^2}{2}f''(x) + O(h^3)\\ f(x+2h) &= f(x) + 2hf'(x) + \frac{4h^2}{2}f''(x) + O(h^3) \end{aligned}$$
After you find the solution $$\alpha_1 = -\frac{3}{2h}, \alpha_2 = \frac{2}{h}, \alpha_3=-\frac{1}{2h}$$
You can find the error between the true and approximated one of each point is of order 2 $$|f'(x)-g(x)| = O(h^2)$$
(Note that it's 2 instead of 3 since the coefficients cancel one $h$)If we choose 2-norm as the error metric, then the error is $$e_k = \left(\sum_{j=1}^n (f'(x_j)-g(x_j))^2\right)^{\frac{1}{2}} = \sqrt{n (O(h^2)^2)} = 2^{\frac{k}{2}}O((2\pi)^2 2^{-2k}),$$ where $n=2^k, k=3 \dots 10$ is the number of mesh points and $h$ is the mesh size $h=\frac{2\pi}{n}$
Then the rate of convergence in the 2-norm case is $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|} = 2^{-\frac{3}{2}},$$ which coincides with the numerical result.
In addition, you can check the case of using mean squared error as the error metric, $$e_k = \frac{1}{n}\left(\sum_{j=1}^n (f'(x_j)-g(x_j))^2\right) = O(h^2)^2 = O((2\pi)^4 2^{-4k}),$$ where $n=2^k, k=3 \dots 10$ is the number of mesh points and $h$ is the mesh size $h=\frac{2\pi}{n}$
Then the rate of convergence in mean squared error is $$\lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|} = 2^{-4}$$
(Use MSE as the error metric)
WolframAlpha is a good tool to calcualte integral like this.
Note that the implementation of Gauss-Legendre requires the parameter $n$ as the number of points (and this $n$ also means the order of used Legendre polynomials).
So it should be "p_roots(n)" instead of "p_roots(n+1)"!
(Ref. Wikipedia)