Problem implementing finite difference method at edges of periodic functions
6 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
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.
HTH
6 Comments
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...
See Also
Categories
Find more on Boundary Conditions in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!