On Mon, Jun 22, 2009 at 10:12 AM, John B. Thoo
You can just think of gain as a constant that multiplies the (Fourier
or Chebyshev) coefficients. If the gain were 1.0 all the way across
you would not be doing any filtering. The filters in the graphs are
low-pass filters so the gain is close to 1.0 for low frequencies and
approaches 0.0 for high frequencies.
You can think of a discontinuity in your solution as a feature with
very high (infinite) frequency. So by shrinking the high frequencies
you are shrinking the influence of the discontinuity, and unfortunately
any other high frequency feature in your solution as well. This is where you need to balance your grid resolution with the order of your filter and the strength of your discontinuity to resolve the features you care about without getting nonsense solutions.
On Jun 19, 2009, at 12:08 PM, Joshua Stults wrote:
I'll second that, you *have* to do something about discontinuities in
your solution, or your solution will become nonsense (blow-up) in
rather few time-steps (depending on your grid resolution compared to
the strength of the discontinuity and any viscosity, artificial or
otherwise).
On Fri, Jun 19, 2009 at 2:23 PM, John B. Thoo<address@hidden> wrote:
Hello, Thomas and everyone.
On May 13, 2009, at 8:53 PM, Thomas Shores wrote:
I'm not going to dissect your code, but I will include a template
script that reduces your problem to making a well-formed definition
of the complex rhs f(z,t). It's at the bottom of this email, and
if you remark/deremark the appropriate lines, you'll get each of
the first two examples you proposed.
However, I have a few comments about your method, which is a
separation of variables type argument with complex exponentials
whose time-dependent coefficients are just the Fourier coefficients
of the solution:
1. There is an immediate problem with truncation m=-M..M of the
double series in the indices k, m which comes originally from a
double series in k and ell, with k+ell=m. The problem is that for
all indices k except k=0, you will be missing some terms with index
k-m, which compounds the mathematical truncation error. It also
means you have to be careful about limits when you calculate the
sums involved in the rhs.
An alternative way to deal with shocks is upwind schemes, which I think
some other folks on the list already pointed you towards, with an upwind scheme there's none of this "fiddling" with filters.
2. There is a question as to why you are using a Fourier method at
all for Burgers equation. This nonlinear hyperbolic problem will
generate shocks in the solution at a minimum time T = -min u_0(x),
where u_0(x) = u(x,0) is the initial condition and the derivative
u_0(x)' is somewhere negative. At this time, the problem becomes
ill-posed and Fourier expansions are not helpful. It is this very
difficulty that is the starting point for the modern theory of
"numerical methods for conservation laws" (see, e.g., LeVeque's
monograph of the same name.)
You can apply a psuedospectral method to this problem, but you need to
filter your solution. For instance checkout the graph at the bottom
of this post:
http://j-stults.blogspot.com/2009/06/chebyshev-spectral-method-shock.html
It's a Chebyshev psuedospectral derivative of a square wave with a
Gaussian blip added (sound familiar?). With this (or any)
psuedospectral method, a discontinuity *anywhere* in your domain will
cause ringing (Gibbs phenomenon) *everywhere* in your domain (because
a spectral method uses global interpolating functions rather than
local ones).
The write-up I linked to above uses a Butterworth filter on the series
coefficients that cleans things up rather nicely, but it still has
some overshoots around the discontinuities, this approach will not
give you the total variation diminishing (TVD) guarantees of the more
appropriate upwind schemes that Thomas mentioned.
Check out the papers by Gotlieb(s) on the equivalence of filtering and
adding a numerical, or artificial viscosity, also about psuedospectral
methods applied to discontinuous problems.
I don't know if this is of any continuing interest, but I want to
bring it (almost) to a close.
Burgers' equation has been of interest for over half a century, and I
think it shows no signs of losing it's appeal. You can learn a lot
from solving Burgers' equation in various ways. It's neat to see how
easy it can be to do it in Octave.
If you are serious about learning spectral methods, I would highly
recommend the book by Trefethen: "Spectral Methods in Matlab", which
will be very relevant since you are using Octave, and it's probably in
your school's library.
Also
Boyd's "Chebyhev and Fourier Spectral Methods".
You can download a copy of Boyd's book so that's nice, my only
criticism would be that the organization is a bit unconventional so it
can be hard to find what you're looking for if you aren't already
pretty familiar with his text.
A more specialized text that does address the issues you run into with discontinuities is volume II of "Spectral methods", by Canuto, Hussaini, Quarteroni and Zang. I'd recommend this one if you are interested in fluid dynamics applications.
After longer than it should have taken, I've finally written a short
code to solve u_t + f'(u)*u_x = 0 using a "pseudospectral" method
avoiding the immediate problem with truncation m=-M..M of the double
series in the indices k, m. (I understand your objections to using a
Fourier method, but that's what I've been directed to do, so I'm
giving it my best.)
Here is the my code. I would appreciate any comments you may have
about it as I'm still very much learning. A few notes:
1) I had to use .' in a couple of places where I hadn't expected,
namely, in the lines
u_x = ifft (ifftshift (i*k.*fftshift (fft (u)).'));
u_t = real (-fprime (u).*u_x.');
Obviously I still don't understand when something is outputted as a
column vector or a row vector. Is there a rule of thumb that can
guide me?
2) The solution does not seem to travel a distance one in time one
when I set f'(u) = 1. Should this be a red flag?
%%------------begin code---------------