numerical methodspartial differential equations
In my notes I have the following definition
If a numerical solution converges to the solution of PDE, then the order of convergence is $$ R = \frac{\log_2 \| \frac{e_{new}}{e_{old}}\|}{\log_2 \| \frac{\Delta x_{new}}{\Delta x_{old}} \|} $$
I'm trying to understand the notation of this formula. What do they mean by $e_{new},\Delta x_{new}, $ etc ? Take for example $v_t+v_x=0$ and suppose we approximate solution with any scheme, for example Crank-Nicolson or leapfrog. How can I compute $R$?
The implementation of two-step methods such as the Leapfrog scheme requires more data vectors than the implementation of one-step methods. One should not forget to update all the data vectors while iterating. Lastly, the initialization must be performed in an upwind fashion. Here is a revised code:
%%%%Here I am finding u_k^1, I used Euler's method to do so.
u1 = u0;
for k=2:N+1
u1(k)=u0(k)-(dt/dx)*(u0(k)-u0(k-1));
end
u1(1)=u1(2);
unm1=u0;
un=u1;
unp1=u1;
for n=1:tsteps
for i=2:N
unp1(i)=unm1(i)-(dt/dx)*(un(i+1)-un(i-1));
end
t=t+dt;
unm1=un;
un=unp1;
plot(x,un,'bo-')
end
with the following output
For the case of $v_t +v_x = 0$, the time step is related to the mesh size via a fixed value of the Courant number $\Delta t/\Delta x$. Therefore, error is analyzed in terms of the mesh size only. As specified in the comments, people like to use powers of two for the number of points $N$. This gives the following type of error plot:
The $x$-axis is the mesh size $\Delta x$ in log scale, the $y$-axis is the error $e$ in log scale. If $e = C (\Delta x)^R$, then we have $\log e = R \log\Delta x + \log C$, i.e. the error plot is a line with slope $R = \frac{\text d \log e}{\text d \log \Delta x}$ in log-log coordinates. For example, one can use the $L^2$-error \begin{aligned} e = \| v_N - v_\text{ref}\| &= \sqrt{\int_0^1 (v_N(x) - v_\text{ref}(x))^2\, \text d x} \\ &\simeq \sqrt{ \Delta x \sum_{j=1}^{N-1} (v_j - v_\text{ref}(x_j))^2 } \end{aligned} but other choices are possible. On the figure, we can see that the slope of the error plot converges towards a fixed value as $\Delta x \to 0$, which is the order of convergence $R \approx 1.57$ measured numerically. To reach the theoretical value of the convergence order, it is important to perform the error measurements by using smooth signals! Then, the estimation is all the more precise as the mesh size is small. A too small mesh size could lead to round-off errors and too expensive computations. Therefore, a compromise is made.
Best Answer
This formula assumes that for the major part, $$e=C\,(Δx)^R.$$ Then one can use the result for two different step sizes to eliminate $C$.
Note that the difference $e$ to the solution is a difference of functions, comparing the exact solution, which is a function, against the implied interpolated function of the discretized solution. This is easiest done within a finite-element framework. In the most simple example you get piecewise linear functions.
Then $\|e\|$ is a function norm. If you reduce it to the samples of the numerical solution and the corresponding samples of the exact solution, you need to compute the norm via a suitable discretization of the function space norm. In the easiest case of the supremums norm you get the maximum norm of the vector without additional factors, in the $L^p$ norms you need to replace the integral with Riemann sums.
With all this then compute