MATLAB Answers

Problem implementing finite difference method at edges of periodic functions

28 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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!