% ks_cnab2.m - solution of Kuramoto-Sivashinsky equation by CNAB2. % u_t = u*u_x - u_xx - u_xxxx, periodic BCs on [0,L] % computation is based on v = fft(u), so linear term is diagonal % compare p27.m in Trefethen, "Spectral Methods in MATLAB", SIAM 2000 % AK Kassam and LN Trefethen, July 2002 % Parameters: space, time, discretization Lx = 32*pi; % length of domain N = 512; % number gridpoints dt = 1/16; % time step T = 100; % integration time x = Lx*(1:N)'/N; % spatial grid % Set initial condition u = cos(x/16).*(1+sin(x/16)).*(1+0.01*sin(x/32)); v = fft(u); % Precompute various CNAB2 quantities: q = (2*pi/Lx)*[0:N/2-1 0 -N/2+1:-1]'; % wave numbers L = q.^2 - q.^4; % Fourier multipliers for linear part of SH eqn D = i*q; % differentiation operator g = -0.5*D; % -1/2 d/dx, for -1/2 d/dx(u^2) = -u u_x % Initialize time-stepping: Nv = g.*fft(u.^2); Nvp = Nv; LR = (1 + (dt/2)*L); LL = (1 - (dt/2)*L); LL_inv = LL.^(-1); nmax = round(T/dt); nplt = floor((T/100)/dt); uu = u; tt = 0; for n = 1:nmax t = n*dt; Nvp = Nv; Nv = g.*fft(real(ifft(v)).^2); v = LL_inv.*(LR.*v + 1.5*dt*Nv - 0.5*dt*Nvp); %v = LL_inv.*(LR.*v); if mod(n,nplt)==0 u = real(ifft(v)); uu = [u, uu]; tt = [tt,t]; end end %min(min(uu)) %max(max(uu)) % Plot results: surf(tt,x,uu), shading interp, lighting phong, axis tight view([90 90]), colormap(jet); set(gca,'zlim',[0 50]) %light('color',[1 1 1],'position',[-1,2,2]) %caxis([-6 2]) %material([0.30 0.60 0.60 40.00 1.00]); %material('dull'); xlbl('t') ylbl('x') %set(gca, 'xtick', []) %set(gca, 'ytick', []) %set(gca,'xcolor',[0.8 0.8 0.8]) %set(gca,'ycolor',[0.8 0.8 0.8]) colorbar