Asked by 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

Answer by Bjorn Gustavsson
on 7 Oct 2019

Edited by 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

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));

end

I get an extremely strange output, looking like this:

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
on 7 Oct 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.