Problem implementing finite difference method at edges of periodic functions

I am working on a problem where I want to use ode45 on the KS equation . To do so, I am employing the method of lines by semidiscretising the equation in the spacial dimension. The code has periodic boundary conditions, However, I am running into problems when trying to test the finite difference expressions for the second and fourth order derivatives.

My code for the derivatives are as follows (using central difference approximation):

u_xx = ( u(wrapN(i+1,N)) - 2*u(i) + u(wrapN(i-1,N)) ) / (h^2);
u_xxxx = ( u(wrapN(i+2,N)) - 4*u(wrapN(i+1,N)) + 6*u(i) - 4*u(wrapN(i-1,N)) + u(wrapN(i-2,N))) / (h^4);

where wrapN is a helper function to help wrap the indices at the boundaries of u to the start/end of the input matrix (and returns the wrapped value fine):

wrapN = @(x, n) (1 + mod(x-1, n));

and N is simply defined as the length of the input array.

However, when I create a test case using a simple cosine function, I get inaccurate results. When I use the input:

x = -pi:0.01:(pi-0.01);
u = cos(x);

Calculating u_xx centered at i=1 gives 1.3692, rather than 1 as expected.

Calculating u_xxxx at i=1 gives an even more incorrect value, -7.8937e+03, and additionally u_xxxx at i=N gives 4.7062e+03 rather than a similar value as expected for a periodic function.

The expressions are otherwise well behaved for values outside those used by the boundary calculations

Answer 1

[Bjorn Gustavsson]

What I'd usualy have to do in this situation is to extract the differentiation matrices - just to check that all elements become corrrect.

Another point that looks problematic is that your x is not giving you a uniform vector over the periodic intervall [-pi pi).

When I try a more uniform vector x:

x = linspace(-pi,pi,321);
x = x(1:end-1);

then at least the second difference of your test comes out right.

HTH

[George Lu]

Hi, for some reason when I use this as my domain, I get periodic behaviour as expected, but my boundary values are 3.85 and -14.86 for u_xx and u_xxxx, respectively. Additionally, the function output is no longer a plot of 0.5*sin(2x), as expected.

When I use the uniform x, my output function looks like this:

When I expect a shape like this (from the influence of a u*u_x term, discretized as 1/6h * (u_i+1 + u_i + u_i-1)(u_i+1 - u_i-1)), which the domain that has a discontinuity at the end gives me:

Any idea what causes this?

[Bjorn Gustavsson]

Not really, but are you sure that the curve is that discontinuous? The gradient at the edges are rather steep, if you wrap these around would it become neatly continuous? What if you try something like:

x_long = [x,x+2*pi]; % Just to cover 2 periods
u_long = [u,u]; % For the final step in time?
subplot(2,1,1),plot(x_long,u_long)
subplot(212),plot(x_long,gradient(u_long,x_long))

Does that give you jumps and discontinuities...

[George Lu]

Not that the curve is discontinuous, rather that the use of a 0.01 point spacing likely caused errors from a mismatch in values at either end of -pi and pi, which should be resolved with linspace as you posted. But for some reason using linspace to make the domain creates a completely different output function, for reasons completely beyond me

If it helps, this is my full calculation of the output:

dudt = zeros(N,1);
    for i = 1:N
        u_xx = ( u(wrapN(i+1,N)) - 2*u(i) + u(wrapN(i-1,N)) ) / (h^2);

        u_xxxx = ( u(wrapN(i+2,N)) - 4*u(wrapN(i+1,N)) + 6*u(i) ...
                     - 4 * u(wrapN(i-1,N)) + u(wrapN(i-2,N))) / (h^4);
        uu_x = 1/6 * ( u(wrapN(i-1,N)) + u(i) + u(wrapN(i+1,N)) ) ...
                        * ( u(wrapN(i+1,N)) - u(wrapN(i-1,N)) ) / h;
        dudt(i) = - u_xx - u_xxxx - uu_x;
        if any([i==1, i==N])
            disp(u_xx)
            disp(u_xxxx)
        end
    end

where i plot the result of dudt afterwards

[George Lu]

Addendum: When I increase the amount of points in the linspace to 100001, the resulting function looks like sin(2x) as expected. However, that is only for a value of h=0.01 as before. When I corrected that and set h to the average of the distance between x points:

h = 0;
for i = 2:N
    h = h+1/(N-1) * (x(i)-x(i-1));
end

I get an extremely strange output, looking like this:

[Bjorn Gustavsson]

Instead of calculating h in a loop, why not do it the matlab way:

h = mean(diff(x));
% Then check that you have uniform spacing:
sigma_h = std(diff(x));

That way you need to think way less about details like loop-variables and where this and that...

[George Lu]

I've modified h to be calculated that way, and the standard deviation is quite small (h is on the order of 10^-5, sigma_h is on the order of 10^-16), but there is no change in my output. Do you mind showing me what your output is?