To calculate and plot L1 norm and Linf norm of my Finite difference approximation as well as find the order of accuracy of the error norms

Hi,

My overall problem statement is below:

On the one-dimensional uniform grids with 𝑁𝑗 points given by 𝑥𝑖 = 𝑖−1 /N𝑗−1 where 1 ≤ 𝑖 ≤ 𝑁𝑗 , 𝑁𝑗 = 2 𝑗 + 1, and 1 ≤ 𝑗 ≤ 13, and using the function 𝜙(𝑥) = cos (2𝜋𝑥), do the following:

(a) For each of the following approximations, present a plot of the logarithm of the two error norms (𝐿1- and 𝐿∞-norms) versus the logarithm of the grid spacing (i.e., log(error norms) vs log(∆x)), and determine the order of accuracy of the error norms:

ii. The second-order accurate central-difference approximation assuming no periodicity, for which first-order accurate one-sided finite-difference approximations are used near the boundaries.

iii. The fourth-order accurate central-difference approximation assuming periodicity. For all approximations, compute the order of accuracy for both the 𝐿1- and 𝐿∞-norm of errors 𝑒𝑖 , defined as 𝑒𝑖 = (1 /2𝜋) | 𝛿𝜙 /d𝑥 |𝑖 − 𝑑𝜙/ 𝑑𝑥 |𝑖|

My Code for ii is below. it gives me the below graph. It does not make since, error should drop as grid spacong decreases.

I expect a graph similar to below:

Come someone please help me see if my code is calculating something incorrectly or summing a variable up in a incorrect direction. I cannot correct this for the life of me.

Any help is greatly appreciated.

clear 
clear all
format long

 for j=1:13
     Nj(j) = ((2.^j)+1)';
     for i = 1:Nj(j)
         xi(j,i) = ((i-1)./(Nj(j)-1));
         delx(j,i) = 1/(2^j)';
     end
 end


f = @(xi) cos(2*pi*xi);
der = -2*pi*sin(2*pi*xi);
y = f(xi);
x = length(xi);
dynoper = zeros(j,x);
count = sum(xi ~= 0, 2);


%Assuming Non Periodic
for l=1:j
    for k=1:x
        switch k
            case 1
                dynoper(l,k) = (y(l,(k+1))-y(l,k))./(delx(l,2));
            case count(1)+1
                dynoper(l,k) = (y(l,k)-y(l,(count(l))))./(delx(l,2));
            otherwise if  k <= (count(l))
                    dynoper(l,k) = (y(l,(k+1))-y(l,(k-1)))./(2.*delx(l,2));
                      else
                    dynoper(l,k) = 0;
            end
        end
    end
end

%NON periodic error
for l=1:j
    for k=1:x
    noperei(l,k) = (dynoper(l,k)-der(l,k))./(2.*pi);
    end
end

%NONperiodic L1 and Linf
L1meannoper = sum(noperei, 2) ./ sum(noperei ~= 0, 2);
graphdx = delx(1:j,1);
Linfnoper = max(noperei,[],2,"omitmissing");


%Plotting the actual derivative and the approximate derivative
figure;
hold on
plot (xi(l,(1:k)), dynoper(l,(1:k)),'b');
plot (xi(l,(1:k)), der(l,(1:k)), 'g');%
grid on
legend ('dphi (Appoximate)', 'dphi (True)')
xlabel ('x spacing')
ylabel ('Y Values at x')
title ('Comparing approximate dphi to actual dphi assuming NO periodicity')

%Plottting the error assuming NO periodicity
figure;
plot (xi(6,(1:k-1)),noperei(6,(1:k-1)));
hold on
plot (xi(l,(1:k-1)),noperei(l,(1:k-1)));
grid on
ylabel('Error Magnitude')
title ('Total Error as grid spacing increases assuming NO periodicity')
legend ('J=6','J=13')

%Plottting the L1 error assuming NO periodicity
figure;
plot (graphdx,L1meannoper, 'r');
grid on
xlabel('DelX')
ylabel('L1Norm (Mean)')
title ('Compare L1 with increasing grid spacing assuming NO periodicity')


%Plottting the log(L1) error assuming NO periodicity
figure;
loglog(graphdx,L1meannoper, 'r');
xlabel('Log(DelX)')
ylabel('Log(L1Norm (Mean))')
title ('Compare Log(L1) with increasing grid spacing assuming NO periodicity')
grid on

%Plottting the log(Linf) error assuming periodicity
figure;
loglog(graphdx,Linfnoper, 'r');
xlabel('Log(DelX)')
ylabel('Log(Linf) (Max))')
title ('Compare Log(Linf) with increasing grid spacing assuming periodicity')
grid on

Answer 1

[Umar]

Hi @Brian,

Your task involves analyzing the accuracy of finite-difference approximations for the derivative of the function ( phi(x) = cos(2\pi x) ) on a one-dimensional uniform grid. Specifically, the goal is to compute and plot the ( L_1 ) and ( L_infinity ) error norms against the grid spacing, while determining the order of accuracy for both norms. The provided MATLAB code by you for the second-order accurate central-difference approximation is yielding unexpected results, prompting a review and correction of the code. So, after analyzing your comments and code, I made modifications to it which appears to address several key issues that may have contributed to the unexpected results in yours original implementation. Below, I will outline the primary improvements made in the modified code and how they contribute to resolving the issues you faced.

Key Improvements in the modified code

Initialization of Variables

The modified code now initializes Nj, delx, L1meannoper, and Linfnoper more clearly, establishing a structured approach to storing grid sizes and error norms. In contrast, your code use of nested loops made it less clear and potentially error-prone when referencing indices across different dimensions.

Grid Point Calculation

The calculation of grid points (xi ) is now more straightforward, directly using vectorized operations rather than nested loops. This minimizes potential indexing errors and enhances readability.

Derivative Approximation

The derivative approximation logic has been simplified. The central difference approximation is explicitly handled with clear conditions for boundary points, which reduces complexity and improves accuracy. Your code had a convoluted structure with unnecessary checks that could lead to indexing errors or incorrect calculations.

Error Calculation

The calculation of normalized errors (noperei) is clearer and uses absolute values directly, ensuring that negative discrepancies do not skew the results. L1 and L∞ norms are computed directly from the absolute values of errors, which enhances clarity and correctness.

Plotting

The plotting section is streamlined for better visualization, making it easier to interpret results. The use of markers in `loglog` plots helps identify data points clearly.

Here are some additional insights that I would like to share with you.

Expected Behavior of Error Norms: As grid spacing (Δx) decreases, one would typically expect both L1 and L∞ norms to decrease, indicating improved accuracy of the derivative approximation. If this behavior is not observed:

• Check for potential issues such as incorrect boundary conditions or incorrect calculations of derivative approximations.

• Make sure that the function being differentiated behaves well (i.e., is smooth) over the grid defined by `xi`.

Order of Accuracy: After obtaining the plots for L1 and L∞ norms against log(Δx), one can determine the order of accuracy by analyzing the slope of the log-log plots. A linear fit can provide insights into how quickly errors decrease as grid resolution increases.

Debugging Tips: If unexpected results persist:

• Print intermediate values (e.g., dynoper, noperei) for different iterations to trace where discrepancies arise.
• Compare outputs from both codes for specific values of `j` to pinpoint where differences occur.

Modified Code

Here is the revised MATLAB code that addresses the issues in the original implementation:

clear; 
clc; 
format long;

% Define the function and its derivative
f = @(x) cos(2*pi*x);
df = @(x) -2*pi*sin(2*pi*x);

% Initialize parameters
max_j = 13;
Nj = zeros(1, max_j);
delx = zeros(1, max_j);
L1meannoper = zeros(1, max_j);
Linfnoper = zeros(1, max_j);

% Loop over j to create grids and compute errors
for j = 1:max_j
  Nj(j) = 2^j + 1; % Number of grid points
  xi = (0:Nj(j)-1) / (Nj(j)-1); % Grid points
  delx(j) = 1 / (2^j); % Grid spacing
  y = f(xi); % Function values
  der = df(xi); % True derivative values

    % Initialize the derivative approximation array
    dynoper = zeros(1, Nj(j));

    % Central difference approximation (non-periodic)
    for k = 1:Nj(j)
        if k == 1
            dynoper(k) = (y(k+1) - y(k)) / delx(j); % Forward difference
        elseif k == Nj(j)
            dynoper(k) = (y(k) - y(k-1)) / delx(j); % Backward difference
        else
            dynoper(k) = (y(k+1) - y(k-1)) / (2 * delx(j)); % Central difference
        end
      end

      % Calculate the error
      noperei = (dynoper - der) / (2 * pi); % Normalized error

      % Compute L1 and Linf norms
      L1meannoper(j) = sum(abs(noperei)) / Nj(j); % L1 norm
      Linfnoper(j) = max(abs(noperei)); % L∞ norm
  end

% Plotting the results
figure;
loglog(delx, L1meannoper, 'r-o');
xlabel('Log(Δx)');
ylabel('Log(L1 Norm)');
title('Log(L1 Norm) vs Log(Δx) assuming NO periodicity');
grid on;

figure;
loglog(delx, Linfnoper, 'b-o');
xlabel('Log(Δx)');
ylabel('Log(L∞ Norm)');
title('Log(L∞ Norm) vs Log(Δx) assuming NO periodicity');
grid on;

% Display the results
disp('Grid Spacing (Δx):');
disp(delx);
disp('L1 Norms:');
disp(L1meannoper);
disp('L∞ Norms:');
disp(Linfnoper);

Please see attached.

Explanation of the Code

Function Definitions: The function ( \phi(x) ) and its derivative are defined using anonymous functions for clarity and ease of use.

Grid Generation: The grid points ( x_i ) are generated based on the formula provided, ensuring that the spacing ( \Delta x ) is calculated correctly.

Finite-Difference Approximations: The code implements a second-order accurate central difference for interior points and first-order differences for the boundaries. This is crucial for accurately capturing the behavior of the derivative near the edges of the domain.

Error Calculation: The error is computed as the difference between the approximate derivative and the true derivative, normalized by ( 2\pi ) to match the problem statement.

Norm Calculation: The ( L_1 ) and ( L_\infty ) norms are calculated for each grid size, allowing for a comparison of accuracy as the grid is refined.

Plotting: The results are plotted on a logarithmic scale to visualize the relationship between the error norms and grid spacing, facilitating the determination of the order of accuracy.

This modified code should yield the expected results, with the error norms decreasing as the grid spacing decreases. The logarithmic plots will help in assessing the order of accuracy of the finite-difference approximations. If you have any further questions or need additional modifications, please feel free to ask.

[Umar]

Hi @Brandon,

Please let me know if you need any further assistance.

[Brian]

Hello @Umar

This was a great answer, thank you so much.

The slope I get for the L1 mean error not=rm is around 1.8, which is roughly 2 which is what to expect for a second order approximation. So you are answer was of great value, thank you.

One follow up question,

If I was to do the same thing as above but assume periodic boundaries, how would you edit the central difference approximation to handle periodic boundaries?

My attempt to do this is below and is giving me a slope of the L1mean at around .8. Which is not correct, it should be around 2.

    % Central difference approximation (Periodic)
    for k = 1:Nj(j)
        if k == 1
            dyper(k) = (y(k+1) - y(Nj(j)-1)) / delx(j); % Periodic start point
        elseif k == Nj(j)
            dyper(k) = (y(k-k+2) - y(k-1)) / delx(j); % Periodic End Point
        else
            dyper(k) = (y(k+1) - y(k-1)) / (2 * delx(j)); % Central difference
        end
      end

[Umar]

Hi @Brian,

You mentioned, “One follow up question,If I was to do the same thing as above but assume periodic boundaries, how would you edit the central difference approximation to handle periodic boundaries? My attempt to do this is below and is giving me a slope of the L1mean at around .8. Which is not correct, it should be around 2.”

I did examine your current implementation, please see my brief analysis summary below.

Boundary Conditions: For k = 1, you correctly reference y(Nj(j)-1) as it wraps around to the last element. However, for k = Nj(j), there seems to be an error in your index calculation: y(k-k+2) is not a valid expression. Here is a revised version of your code that properly handles periodic boundaries:

% Central difference approximation (Periodic)
for k = 1:Nj(j)
  if k == 1
      dyper(k) = (y(2) - y(Nj(j)))/(2 * delx(j)); % Corrected for         periodicity at start
  elseif k == Nj(j)
      dyper(k) = (y(1) - y(Nj(j)-1))/(2 * delx(j)); % Corrected         for periodicity at end
  else
      dyper(k) = (y(k+1) - y(k-1)) / (2 * delx(j)); % Central         difference for interior points
  end
end

So, in the above modified snippet code, you will see the following changes made to the code.

Start Point (k == 1): The derivative at the first point should consider the second point (y(2)) and wrap around from the last point (y(Nj(j)))

End Point (k == Nj(j)): The derivative at the last point now correctly references y(1) and wraps around from the second-to-last point (y(Nj(j)-1)).

Now, if you are getting a slope of approximately 0.8 instead of the expected 2, make sure that your values in y are set up correctly and represent the expected function over your domain. Also, implement checks to make sure that delx(j) is not zero to avoid division errors and consider plotting your results or using debugging statements to examine values of dyper during iterations to better understand discrepancies in expected versus actual outcomes.

By carefully revising how you handle boundary conditions in your central difference approximation, you can achieve more accurate results consistent with periodic behavior. Implement these changes, and then verify your outputs against known solutions or analytical results to further validate your approach. If you continue to face issues, reviewing data input or considering numerical stability may also be beneficial.

Hope this helps. Please let me know if you have any further questions.