MATLAB Answers

Problem implementing finite difference method at edges of periodic functions

17 views (last 30 days)
George Lu
George Lu on 7 Oct 2019
Commented: George Lu on 7 Oct 2019
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


Sign in to comment.

Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 7 Oct 2019
Edited: Bjorn Gustavsson on 7 Oct 2019
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.


Show 3 older comments
George Lu
George Lu on 7 Oct 2019
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));
I get an extremely strange output, looking like this:
Bjorn Gustavsson
Bjorn Gustavsson on 7 Oct 2019
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
George Lu on 7 Oct 2019
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?

Sign in to comment.

Sign in to answer this question.