Problem: Stiffness and Spectral Methods



Table of contents:

  (1) Introduction: solitary waves.

  (2) Description of a pseudo-spectral method.

  (3) Initial value and interval to use.

  (4) Code to use/write.

  (5) About ode solvers and other programing hints.

  (6) What to do.



(1) INTRODUCTION: SOLITARY WAVES. =================================



Consider the Korteweg deVries (KdV) equation



              u_t + c*u_x + (0.5*u^2)_x + u_xxx = 0,         (1.1)



where c is a constant. This equation can be used to describe weakly

nonlinear waves in shallow channels.    The equation has the set of

special solutions given by



           u =  3*(p^2)*(sech(0.5*p*z))^2

             = 12*(p^2)*(exp(0.5*p*z) - exp(0.5*p*z))^(-2),  (1.2)



where z = x - s*t - x0, s = c + p^2 and

--- p > 0 is a constant,

--- x0   is a constant,



These solutions are called "solitary waves" because they consist of

a single isolated "bump" in u, moving at speed s. We notice that:

--- The solution maximum occurs at z = 0, where u = 3*p^2 = A.

--- The solution decays exponentially,  as z grows to +/- infinity,

    with the bump having a width of order W = 4/p. In fact:

    For |z| =   W/2, ...... u/A = 0.42 ...

        |z| =   W, ........ u/A = 0.07 ...

        |z| = 2*W, ........ u/A = 0.001 ...

        |z| = 3*W, ........ u/A = 0.000025 ...

        ......................................



Thus,  within a few widths W of z = 0,  this solution (essentially)

vanishes.  Because of this,  if we take two solutions of this type,

and add them, then:

   As long as the two bumps remain far from each other,  their

   sum will also be a solution [this because, in this case the

   nonlinear term has no effect, since at any point one or the

   other solution vanishes].



(2) DESCRIPTION OF A PSEUDO-SPECTRAL METHOD. ======================



Write the equation in the form



         u_t = - (c*u_x + (0.5*u^2)_x + u_xxx),              (2.1)



and now consider solutions  (with periodic boundary conditions)  in

an interval 0 <= x <= Le. Next introduce a uniform mesh x_n, with N

nodes (ASSUME N = 2*M is EVEN), given by the (MatLab) array



                   x = Le*(0:N-1)'/(N-1).                    (2.2)



Then consider the array of values for the solution  u = u(x, t)  at

the nodes x_n (u_n = u_n(t) = u(x_n, t)). This gives the array of N

values u = [u_1; u_2; .... u_N].



In terms of the N-array u, if we use fft's to evaluate derivatives,

equation (2.1) becomes a system of N ode's for the values u_n. This

system of ode's can then be solved any standard ode solver.



  Examples:  u_x ........... becomes real(ifft(k.*fft(u));

             (0.5*u^2)_x ... becomes real(ifft(0.5*k.*fft(u.^2)));



where k is the array  k = i*(2*pi/Le)*[(0:M), (1-M):(-1)]';  (2.3)



The above describes a "pseudo-spectral" method, that can be used to

solve the KdV equation.



(3) INITIAL VALUE AND INTERVAL TO USE. ============================



Consider the solution, with periodic boundary conditions in the in-

terval 0 <= x <= Le, with Le = 40, and the initial condition



               u0 = 3*(sech(0.5*(x - x0))).^2,               (3.1)



where x0 = 2*Le/5.   This corresponds to a solitary wave [see (1.2)

above] of

        amplitude A = 3, width W = 4, and speed s = c+1,     (3.2)



initially located [center of the bump] at x = x0.  Thus, with these

initial conditions, we know exactly what the exact solution is.  We

note that, for c = - 1, the solution does not change in time.



Note: You can change Le and x0. The only real condition is that you

      should have x0 >> W,  and Le-x0 >> W,  so that u0 vanishes at

      both ends of the interval ---- thus u0 satisfies the periodic

      boundary conditions. 



(4) CODE TO USE/WRITE. =============================================



Write a script to solve (1.1) --- with the initial conditions in (3)

--- using the ideas in (2).    Make it so that you can easily change

the values of N = 2*M, c, and the ode solver used [see (5) below for

info about the ode solver  that you will need to implement the stuff

below]. Set it up so that:



a---

The solution is computed for 0 < t < tf, for some tf > 0 that should

be easy to change.

                  A GOOD CHOICE IS tf = 3.



b---

Do the computation in Nt jumps of length ............... dt = tf/Nt.

Again, make Nt easy to change. A GOOD CHOICE IS ........ Nt = 6.

Thus, you will end up computing the solution

for the times



   t_0 = 0;  t_1 = dt;  t_2 = 2*dt; ... t_n = n*dt; ... t_Nt = tf.



c---

After each computation [i.e.:  after you finish with the computation

for time t=t_n]:

   Record the time taken by the ode solver to get from

   t = t_{n-1} to t = t_n, and display the answer.



   Record the number of time steps taken by the ode solver

   to get from t = t_{n-1} to t = t_n, and display the answer.



   Plot the solution  [use "hold on" to make sure the solution

   at prior steps is not erased; also use different colors for

   each time ... and remember to plot the initial value before

   you even start].



(5) ABOUT ODE SOLVERS and OTHER PROGRAMING HINTS. ==================



a---FUNCTION SCRIPT.

You will need a script that the ode solver will call to evaluate the

right hand side in (2.1). For example, create a file called kdvfun.m

and write in it:



% ------- BEGIN kdvfun.m FILE.

function u_dot = kdvfun(t, u, flag, k, c)

if isempty(flag)

   u_dot = - real(ifft((c*k + k.^3).*fft(u) + 0.5*k.*fft(u.^2)));

else

   u_dot = [];

end

% ------- END OF kdvfun.m FILE.



Note that k is the array defined in (2.3) above.  For this script it

is an input.    The function scripts for ode solvers in MatLab allow

extra inputs in addition to t and the function value. The input flag

is the mark beyond where the extra inputs are.



b---ODE SOLVER; etc.

Here follows the core of a script doing the stuff required. Read the

comments to understand how it works.

%

% ==================================================================

% This section contains the parameters that can be changed.

%

c   =    0;% ------- Constant in the equation.

M   =   16;% ------- N = 2*M = number of points in grid.

Le  =   40;% ------- Length of computational interval [0, Le].

tf  =    3;% ------- Final time for computation.

Nt  =    6;% ------- Number of times to plot during calculation.

tol = 1e-4;% ------- Tolerance for ode solver.

COL = {'c', 'r', 'g', 'b', 'm', 'k', 'y'};

% --- COL contains the colors used for the plots.

% --- COL must have Nt+1 colors.

%

% ==================================================================

% Below the grid and initial values are computed.  Note that the ode

% solvers require COLUMN vector inputs, but output rows. The various

% transposes (') that show all over the code it to take care of this

% detail.

%

N  = 2*M;% ---------------------------------- Number of points used.

dt = tf/Nt;% -------------------------------- Time between plots.

x  = Le*(0:N-1)'/N;% ------------------------ Space grid.

u0 = 3*(sech(0.5*(x - (2/5)*Le))).^2;% ------ Initial values.

k  = i*(2*pi/Le)*[(0:M), ((1-M):(-1))]';% --- k Fourier array.

umin = -0.25;% ----------------- Minimum value to use in plot range.

umax =  3.25;% ----------------- Maximum value to use in plot range.

%

% ==================================================================

% Plot the initial conditions and setup figure for further plots.

% Get also some text to show on the command window, so you know all

% is well.

%

tp = 0;% -------------- Will contain current time in calculation.

xp = [ x', Le];% ------ For plotting it is convenient to add the

up = [u0', u0(1)];% --- right point of the interval [0, Le].

%

fprintf(['\n Begin integration with ', num2str(N), ' grid points.'])

fprintf(['\n Current time = ', num2str(tp)])

fprintf('\n Plotting initial conditions.')

figure;% ---- Creates a figure.

plot(xp, up, '-', 'LineWidth', 2, 'Color', COL{1})

xlabel('{\bf x}', 'FontSize', 14)

ylabel('{\bf u}', 'FontSize', 14)

title(['{\bf KdV stiffness demo: N = ', num2str(N), ...

       ', c = ', num2str(c), '}'], 'FontSize', 14)

axis([0 Le umin umax])

grid on;% --- This will put a grid on the plot.

zoom on;% --- This allows zooming with the mouse.

hold on;% --- So new plots are added without erasing prior ones.

pause(1);% -- Pauses execution, so MatLab can plot.

%

% ==================================================================

% Setup the ode solver.

% Note how the odeset command allows control over the tolerance.

% Note how the input k, needed by kdvfun, is given to the solver.

% Here we use ode45, but any other solver would just require a name

% change below.

% The commands tic and tac are used to "time" the solver.

%

options = odeset('RelTol', tol, 'AbsTol', tol/1000);

Y0      = u0;

for n=1:Nt

   tic;

   [T, Y] = ode45('kdvfun', [0, dt], Y0, options, k, c);

   ts = toc;

   Ns = length(T)-1;

   tp     = tp+dt;

   Y0     = (Y(Ns+1, :))';

   up     = [Y0', Y0(1)];

   fprintf(['\n\n Integration number ', num2str(n), ' finished.'])

   fprintf(['\n Current time ..................... ', num2str(tp)])

   fprintf(['\n Number of integration steps ...... ', num2str(Ns)])

   fprintf(['\n Time taken to integrate (secs) ... ', num2str(ts)])

   fprintf('\n Press any key to plot u at current time.')

   pause

   plot(xp, up, '-', 'LineWidth', 2, 'Color', COL{n+1})

   pause(0.5)

end

fprintf('\n\n THE END \n')

% ------------------------ END OF SCRIPT



(6) WHAT TO DO. ====================================================



Use the code you just wrote to explore what happens when you solve

the equations with various values of N, and various solvers.  Try,

for example:



ode45  and N = 32, 64, 128, ...

ode15s and N = 32, 64, 128, ...



What happens as N increases?

What happens as the integration goes forward

     [i.e. from t=0, to t=dt, to t=2*dt ...].

Zoom into the solution [near the extremes of the computing interval,

     where it is supposed to be zero]. What do you see?

     [Note: the code above enables the mosue zoom feature].

What happens when different solvers are used?



Write a short report with your conclusions.

%

%% EOF