28 views (last 30 days)

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

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

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

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.