Answer Happy

Hi, I posted this yesterday and saw know that there was a slight
type error in the section5.m code. I have fixed it and will post
the question again!
(I'm using matlab) - Numerical approximation of solitary
waves.
Run the code kdvdis.m, and study the appearance of solitary
waves from Gaussian initial data for the KdV equation. Use the code
Section5.m to generate approximate solitary waves for the Whitham
equation (see Figure 1). Using a numerical cleaning technique (i.e.
manually zeroing out undesirable parts of the solution), identify
two approximate solitary waves, one of small amplitude, and one of
large amplitude.
kdvdis.m:
%kdvdis.m: Disintegration of an initial Gaussian pulse
%into solutions: u_t + disp*u_{xxx} + uu_x = 0
clear;
N=512;
L = 30;
a = L/(2*pi);
T = 50;
h = 2*pi/N;
dt = .005;
y = a*(h:h:2*pi)';
u = exp(-2*(y-.25*L).^2)/a;
k = [(0:N/2)';(1-N/2:-1)'];
disp = 0.01;
m = disp*a^(-3)*.5*dt*k.^3;
d1 = (1+i*m)./(1-i*m);
d2 = -.5*i*dt*k./(1-i*m);
d3 = .5*d2;
sol = plot(y, a*u, 'Erasemode', 'background');
axis([ 0 L -.3 2]);
zoom;
title('KdV: Disintegration of a Gaussian pulse');
drawnow;
t=0;
while t<T
v = real(ifft(d1.*fft(u) + d3.*fft(u.^2)));
w = real(ifft(d1.*fft(u) + d2.*fft(u.^2)));
for n = 1:3
w = v + real(ifft(d3.*fft(w.^2)));
end
u = w;
t = t + dt;
set(sol, 'ydata', a*u);
end
section5.m:
%isospectral integration of u_t + uu_x + K*u = 0
% K = fourier transform of dis(k)
N = 512; %input('N=')
L = 8;
a = L/(2*pi);
T = 1.1; %input('T=')
h = 2*pi/N;
dt = .005;
y = a*(h:h:2*pi)';
c1 = 2;
u = exp(-c1*a*(y-pi).^2)/a;
k = [(1:N/2)'; (1-N/2:-1)'];
dis = [1; sqrt(tanh(k)./k)];
Dv = [0; k];
m = a^(-3)*.5*dt*dis.^3;
d1 = (1+i*m)./(1-i*m);
d2 = -.5*i*dt*Dv./(1-i*m);
d3 = .5*d2;
sol = plot(y, a*u, 'r', 'Erasemode', 'background');
box off;
axis([ 0 L -.2 1.75]);
title('Whitham Equation, weak dispersion');
text(1, 1.6, 'dispersion relation = (tanh(k)/k)^{1/2}')
text(1, 1.4, 'u_t + uu_x + K*u = 0, u(x,0) = 2\pi
exp(-2(x-4)^2)/8')
t=0;
while t<T
fftu = fft(u);
fftuu = fft(u.^2);
v = real(ifft(d1.*fftu + 0.5*d2.*fftuu));
w = real(ifft(d1.*fftu + d2.*fftuu));
for n = 1:3
w = v + real(ifft(d3.*fft(w.^2)));
end
u = w;
t = t + dt;
set(sol, 'ydata', a*u);
drawnow;
end
0.5 0.2 Time: 40000 0.45 0.1 0.4 0 0,35 0.2 Time: 10000 0.3 0.1 0.25 0,2 0 Stokes number: S-4.46. 0.15 Whitham regime parameters: 0.2 -Time: 0.05 K=1.35;=0.51. 0.1 Whitham number: W=3.78. 0.1 0.05 0 0 -4000 -3000 -2000 -1000 1000 8 10 12 14 16 18 20 Wavelength [/hol x [x/ho] Fig. 1. Left panel: Formation of solitary waves of the Whitham equation from Gaussian initial data. Right panel: Curve fit for the Whitham regime and for the Boussinesq regime to amplitude/wavelength data from Whitham solitary waves. The wavelength is defined as I = n(x)dx. Amplitude [a/ho] o Data Boussinesq regime Whitham regime