help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Trying to solve du/dt = i*u


From: Joshua Stults
Subject: Re: Trying to solve du/dt = i*u
Date: Fri, 19 Jun 2009 15:08:08 -0400

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

> 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---------------
>
> % This script solves the pde  u_t + f'(u)*u_x = 0  using a
> "pseudospectral"
> % method.
>
>
> % set time and space parameters
> t0 = 0;   % initial time
> tf = 1;   % final time
> L = -1;   % left end point
> R = 1;   % right end point
> M = 10;   % number of time intervals
> N = 100;   % number of space intervals
>
> dt = (tf - t0)/M;   % time mesh
> t = t0:dt:tf;   % time interval
> dx = (R - L)/N;   % space mesh
> x = a:dx:b;   % space interval
>
>
> % define initial condition  u0 = u(x,0)
> %u0 = (x <= 0);   % step function: 1 on L, 0 on R
> u0 = exp (-25*x.^2);   % Gaussian hump
>
> global lenu0;
> lenu0 = length (u0);
>
>
> % define the equation
> function speed = fprime (u)
>   speed = u;   % set  fprime(u)  in  u_t + fprime(u)*u_x = 0
> endfunction
>
> function u_t = F (u,t)
> global lenu0;
>
>   speed = fprime (u);
>
>   c = floor (lenu0/2);
>   k = -c:c;
>
>   u_x = ifft (ifftshift (i*k.*fftshift (fft (u)).'));
>
>   u_t = real (-fprime (u).*u_x.');
> endfunction
>
>
> % solve the equation
> u = lsode("F", u0, t).';
>
>
> % plot the solution
> clf
>
> %% static plot
> %plot(x, u(:,1), 'r', x, u(:,2), 'g', x, u(:,3), 'b')
>
> %% animation
> wait = max (1/(tf - t0 + 1), 0.05);   % pause time between frames
>
> for j = 1:M+1
>   plot (x, u(:, j), '-o');
>   title (strcat ('time =', num2str (t(j))));
> %  axis ([-1, 1, -1, 1]);
>   pause(wait);
> endfor
>
> %%------------end code---------------
>
> Thanks for your help.
>
> ---John.
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>

-- 
Joshua Stults
Website: http://j-stults.blogspot.com



reply via email to

[Prev in Thread] Current Thread [Next in Thread]