% % code7b.m -- spectral wave w/ advective nonlinearity % % u + c u u = d u % t x xxxx % % spectral, integrating-factor formulation: U = F[u]; r=4 % % (U exp(-d (ik)^r t)) = exp(-d (ik)^r t) (-0.5*c) (F[u^2]) % t x % %- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - global ik choice = 4; switch choice case {4} c = 1.0; d = -0.00064; r = 4; par = [c,d,r]; L = 2*pi; N = 128; dt = 0.01; Nsteps = 50; Ncycle = 8; blurb = '\bf 4^{th}-order hyper-diffusive regularization'; case {3} c = 1.0; d = -0.00008; r = 4; par = [c,d,r]; L = 2*pi; N = 256; dt = 0.005; Nsteps = 100; Ncycle = 8; blurb = '\bf 4^{th}-order hyper-diffusive regularization'; case {1} c = 1.0; d = -0.00001; r = 4; par = [c,d,r]; L = 2*pi; N = 512; dt = 0.0025; Nsteps = 200; Ncycle = 8; blurb = '\bf 4^{th}-order hyper-diffusive regularization'; end % set grids: wavenumbers & coefficients ik = i*(2*pi/L)*(-N/2:N/2-1)'; fact = exp(-d*dt*(ik.^r)); dx = L/N; xg = (0:dx:L-dx)'; x = [xg ; L]; % IV initialization & FFT switch choice case {5} IVs = exp(-2*(xg-pi).^2); pmin = -1.1; pmax = 1.1; case {1,2,3,4} IVs = sin(xg); pmin = -1.5; pmax = 1.5; end t = 0; ug = IVs; u = [ug ; ug(1)]; %clf; hold on; axis([0 L pmin pmax]); hold on; axis([0 L pmin pmax]); plot(x,u,'r'); plot(x,u,'k.','markersize',7); vg = fftshift(fft(ug))/N; % timestep loop for k = 1:Ncycle for j = 1:Nsteps [vg,tj] = code7a_rk(vg,0,dt,par); vg = vg./fact; t = t + dt; end ug = real(ifft(fftshift(vg))*N); plot(x,[ug ; ug(1)],'b'); plot(x,[ug ; ug(1)]); pause(0.1) end plot(x,[ug ; ug(1)],'r'); plot(x,[ug ; ug(1)],'k.','markersize',7); title(blurb); xlabel('\bf x-axis'); ylabel('\bf u-axis'); hold off deriv = (ug(N/2+2)-ug(N/2))/(2/N)