%% We solve here the coupled dispersive/dissipative system: %% i u_t + u_{xx} + gamma u = n u + f %% i n_t - |D| n = |D| ( |u|^2 ) %% The code will be spectral in nature due to the non-local operator. % Input parameters for numerical integration %gamma = 1.0; A = 1.905; T = 300.0; gammas = .225; gammad = .225; eta = 1.0; N = 32; %Make this input even. i = sqrt(-1); ht = 1e-4; Nt = round(T/ht); % Setting up the period 2*pi spatial domain. dx = pi/N; x = [dx:dx:2*pi]; N1 = size(x,2); xi = [0:N1/2-1 0 -N1/2+1:-1]'; absXI = abs(xi); absXIsquared = (xi).^2; gammavecs = gammas+0*xi; gammavecd = gammad+0*xi; %mchk = ones(N1,1); %mchk(1) = 0; %mchk(N1/2) = 0; % Initial Data and Potential u0 = A*sin(x)'; n0 = 0*sin(x)'; v0 = (fft(u0)); m0 = (fft(n0)); v0flag = 0*v0; m0flag = 0*m0; v0flag1 = 0*v0; m0flag1 = 0*m0; v0flagmin = zeros(N1,50); m0flagmin = zeros(N1,50); v0flagmax = zeros(N1,50); m0flagmax = zeros(N1,50); forces = eta*fft(sin(x)'); forced = 0*fft(sin(x)'); t = zeros(Nt+1,1); flag = zeros(Nt+1,1); flag2 = zeros(Nt+1,1); t(1) = 0; zsch = zeros(Nt+1,1); zdir = zeros(Nt+1,1); mode0sch = zeros(Nt+1,1); mode1sch = zeros(Nt+1,1); mode2sch = zeros(Nt+1,1); mode0dir = zeros(Nt+1,1); mode1dir = zeros(Nt+1,1); mode2dir = zeros(Nt+1,1); zsch(1) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zdir(1) = (dx*sum(abs(ifft(m0)).^2))^(1/2); mode0sch(1) = v0(N1/2); mode1sch(1) = v0(2); mode2sch(1) = v0(N1); mode0dir(1) = m0(N1/2); mode1dir(1) = m0(2); mode2dir(1) = m0(N1); zavesch = 0; zavedir = 0; zflag = 0; flagmax1 = 1; flagmin1 = 1; foomin = 0; Imin = 0; Imax = 0; zz = 2; Tperiod = 0; %Nonlinear potentials and forcing % Solution Operators for Schrodinger and Dirac Es = exp(-i*ht*(absXIsquared - i*gammavecs)); Es2 = exp(-i*ht*(absXIsquared - i*gammavecs)/2.0); Ed = exp(-i*ht*(absXI-i*gammavecd)); Ed2 = exp(-i*ht*(absXI-i*gammavecd)/2.0); % Setting up the numerical method from T. Potter/Trefethen-Kassam % which discretizes the forward linear operator using an Explicit/Implicit % Runge-Kutta scheme. M = 64; Ls = -i*absXIsquared - gammavecs; Ld = -i*absXI - gammavecd; r = exp(2*i*pi*((1:M)-.5)/M); LRs = ht*Ls(:,ones(M,1))+r(ones(N1,1),:); LRd = ht*Ld(:,ones(M,1))+r(ones(N1,1),:); Qs = ht*mean( (exp(LRs/2)-1)./LRs, 2 ); Qd = ht*mean( (exp(LRd/2)-1)./LRd, 2 ); f1s = ht*mean( (-4-LRs+exp(LRs).*(4-3*LRs+LRs.^2)) ./ (LRs.^3), 2 ); f1d = ht*mean( (-4-LRd+exp(LRd).*(4-3*LRd+LRd.^2)) ./ (LRd.^3), 2 ); f2s = ht*mean( (2+LRs+exp(LRs).*(-2+LRs)) ./ (LRs.^3), 2 ); f2d = ht*mean( (2+LRd+exp(LRd).*(-2+LRd)) ./ (LRd.^3), 2 ); f3s = ht*mean( (-4-3*LRs-LRs.^2+exp(LRs).*(4-LRs)) ./ (LRs.^3), 2 ); f3d = ht*mean( (-4-3*LRd-LRd.^2+exp(LRd).*(4-LRd)) ./ (LRd.^3), 2 ); % Iterative Loop for Convergence Checks of Code with $u(x,t) = n(x,t) = e^{it} sin(x)$. % tic % for j = 1:Nt % t = j*ht; % Nvs = RHSs(v0,m0,forces+0*fft(exp(i*t)*(i*gamma*sin(x)'-2*sin(x)'))); % Nvd = RHSd(v0,absXI,forced+0*fft(exp(i*t)*(i*gamma*sin(x)'-sin(x)'-(ifft(absXI.*fft(sin(x)')))))); % as = Es2.*v0+Qs.*Nvs; % ad = Ed2.*m0+Qd.*Nvd; % Nas = RHSs(as,ad,forces+0*fft(exp(i*t)*(i*gamma*sin(x)'-2*sin(x)'))); % Nad = RHSd(as,absXI,forced+0*fft(exp(i*t)*(i*gamma*sin(x)'-sin(x)'-(ifft(absXI.*fft(sin(x)')))))); % bs = Es2.*v0 + Qs.*Nas; % bd = Ed2.*m0 + Qd.*Nad; % Nbs = RHSs(bs,bd,forces+0*fft(exp(i*t)*(i*gamma*sin(x)'-2*sin(x)'))); % Nbd = RHSd(bs,absXI,forced+0*fft(exp(i*t)*(i*gamma*sin(x)'-sin(x)'-(ifft(absXI.*fft(sin(x)')))))); % cs = Es2.*as+Qs.*(2.0*Nbs-Nvs); % cd = Ed2.*ad+Qd.*(2.0*Nbd-Nvd); % Ncs = RHSs(cs,cd,forces+0*fft(exp(i*t)*(i*gamma*sin(x)'-2*sin(x)'))); % Ncd = RHSd(cs,absXI,forced+0*fft(exp(i*t)*(i*gamma*sin(x)'-sin(x)'-(ifft(absXI.*fft(sin(x)')))))); % v0 = Es.*v0 + (Nvs.*f1s + 2*(Nas+Nbs).*f2s+Ncs.*f3s); % m0 = Ed.*m0 + (Nvd.*f1d + 2*(Nad+Nbd).*f2d+Ncd.*f3d); % end % toc % figure; hold on % Evaluation of dispersive-dissipative model Zakharov system: tic for j = 1:Nt t(j+1) = j*ht; Nvs = RHSs(v0,m0,forces); Nvd = RHSd(v0,absXI,forced); as = Es2.*v0+Qs.*Nvs; ad = Ed2.*m0+Qd.*Nvd; Nas = RHSs(as,ad,forces); Nad = RHSd(as,absXI,forced); bs = Es2.*v0 + Qs.*Nas; bd = Ed2.*m0 + Qd.*Nad; Nbs = RHSs(bs,bd,forces); Nbd = RHSd(bs,absXI,forced); cs = Es2.*as+Qs.*(2.0*Nbs-Nvs); cd = Ed2.*ad+Qd.*(2.0*Nbd-Nvd); Ncs = RHSs(cs,cd,forces); Ncd = RHSd(cs,absXI,forced); v0 = Es.*v0 + (Nvs.*f1s + 2.0*(Nas+Nbs).*f2s + Ncs.*f3s); m0 = Ed.*m0 + (Nvd.*f1d + 2.0*(Nad+Nbd).*f2d + Ncd.*f3d); zsch(j+1) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zdir(j+1) = (dx*sum(abs(ifft(m0)).^2))^(1/2); flag(j) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2) ... + (dx*sum(abs(ifft(m0)).^2))^(1/2) ... - (dx*sum(abs(ifft(v0flag)).^2+ abs(ifft(xi.*v0flag)).^2))^(1/2) ... - (dx*sum(abs(ifft(m0flag)).^2))^(1/2); %flag2(j) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2) ... %+ (dx*sum(abs(ifft(m0)).^2))^(1/2); flag2(j) = (zsch(j+1)^2 + zdir(j+1)^2)^(1/2); mode0sch(j+1) = v0(N1/2); mode1sch(j+1) = v0(2); mode2sch(j+1) = v0(N1); mode0dir(j+1) = m0(N1/2); mode1dir(j+1) = m0(2); mode2dir(j+1) = m0(N1); if mod(j,1e1) == 0 flagconv = abs((dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2)... ... + (dx*sum(abs(ifft(m0)).^2))^(1/2) ... - (dx*sum(abs(ifft(v0flag1)).^2+ abs(ifft(xi.*v0flag1)).^2))^(1/2) ... - (dx*sum(abs(ifft(m0flag1)).^2))^(1/2)); if flagconv < 1e-10 zflag = 1.0; zavesch = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zavedir = (dx*sum(abs(ifft(m0)).^2))^(1/2); zz = 0; fprintf('Converged Orbit') break end v0flag1 = v0; m0flag1 = m0; end v0flag = v0; m0flag = m0; if j > 10 if (flag(j)>=0 && flag(j-1)<0) flagmin1 = flagmin1 + 1; tflagmin(flagmin1+1) = j; v0flagmin(:,flagmin1+1) = ifft(v0); %/(max(abs(v0))); m0flagmin(:,flagmin1+1) = ifft(m0); %/(max(abs(m0))); zperschmin(flagmin1+1) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zperdirmin(flagmin1+1) = (dx*sum(abs(ifft(m0)).^2))^(1/2); for k = 1:flagmin1 foomin1(k) = max(abs(v0flagmin(:,k)-v0flagmin(:,flagmin1+1))); barmin1(k) = max(abs(m0flagmin(:,k)-m0flagmin(:,flagmin1+1))); end [foomin,Imin] = min(foomin1(1:flagmin1)+barmin1(1:flagmin1)); %plot(x,real(ifft(v0)),'--b',x,real(ifft(m0)),'-b',x,real((v0flagmin(:,Imin))),'--r',x,real((m0flagmin(:,Imin))),'-r'); %ylim([-10 10]); %F = getframe; end if (flag(j)<=0 && flag(j-1)>0) flagmax1 = flagmax1 + 1; tflagmax(flagmax1+1) = j; v0flagmax(:,flagmax1+1) = ifft(v0); %/(max(abs(v0))); m0flagmax(:,flagmax1+1) = ifft(m0); %/(max(abs(m0))); zperschmax(flagmax1+1) = (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zperdirmax(flagmax1+1) = (dx*sum(abs(ifft(m0)).^2))^(1/2); for k = 1:flagmax1 foomax1(k) = max(abs(v0flagmax(:,k)-v0flagmax(:,flagmax1+1))); barmax1(k) = max(abs(m0flagmax(:,k)-m0flagmax(:,flagmax1+1))); end [foomax,Imax] = min(foomax1(1:flagmax1)+barmax1(1:flagmax1)); % plot([-N1/2+1:1:N1/2],real(ifftshift(v0)),'--b',[-N1/2+1:1:N1/2],real(ifftshift(m0)),'-b',[1:1:N1],real(ifftshift(v0flagmax(:,Imax))),'--r',[1:1:N1],real(ifftshift(m0flagmax(:,Imax))),'-r'); % ylim([-25 25]); % F = getframe; end if (flagmin1 > 2 && flagmax1 > 2) if (foomin < 5e-3 && foomax < 5e-3) Tperiod = tflagmin(flagmin1+1) - tflagmin(Imin); zavesch = (mean(zperschmax(Imax:flagmax1+1)) + mean(zperschmin(Imin:flagmin1+1)))/2.0; zavedir = (mean(zperdirmax(Imax:flagmax1+1)) + mean(zperdirmin(Imin:flagmin1+1)))/2.0; zflag = 1; zz = 1; fprintf('Periodic Orbit') break end end end %if mod(j,1e3) == 0 % plot([1:1:N1],real(ifftshift(v0)),'--b',[1:1:N1],real(ifftshift(m0)),'-b'); % ylim([-25 25]); % F = getframe; %end if t(j) > T/5 %if mod(j,1e0) == 0 zflag = zflag + 1; zavesch = zavesch + (dx*sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); zavedir = zavedir + (dx*sum(abs(ifft(m0)).^2))^(1/2); %end end % if t(j) > 3*T/5 % if mod(j,1e3) == 0 % plot(zdir(j),zsch(j),'.b','LineWidth',2) % legend('Dirac vs. Schrödinger'); % set(gca,'fontsize',16); % xlabel('Dirac','FontSize',16); % ylabel('Schrödinger','FontSize',16) % xlim([0 20.0]); % ylim([0 15.0]); % F = getframe; % end % end % if mod(j,1e3) == 0 % plot(x,real(ifft(v0)),'--b',x,real(ifft(m0)),'-b','LineWidth',2); % ylim([-10.0 10.0]); % legend('Schrodinger','Dirac'); % set(gca,'fontsize',16); % xlabel('$x$','interpreter','latex','FontSize',16); % ylabel('Computed Solution (Real Part)','FontSize',16) % F = getframe; % end mm0 = ifftshift(m0); if (abs(mm0(1)) > 1 && abs(mm0(N1)) > 1) fprintf('Edge Fourier Modes Too Large, Computation Terminated') break end if mod(j,1e3) == 0 ht*j end end toc hold off figure; plot(t(2:j),flag2(2:j),'--r',t(2:j),zsch(2:j),'-ob',t(2:j),zdir(2:j),'-xg','LineWidth',2) legend('Total Energy','Schrodinger Energy','Dirac Energy'); set(gca,'fontsize',16); xlabel('$t$','interpreter','latex','FontSize',16); ylabel('Energy','FontSize',16) t4 = t(2:j); te4 = flag2(2:j); figure; plot(zdir(round(2*j/4):j),zsch(round(2*j/4):j),'-b',zdir(j),zsch(j),'xk','LineWidth',2) legend('Dirac vs. Schrödinger'); set(gca,'fontsize',16); xlabel('Dirac','FontSize',16); ylabel('Schrödinger','FontSize',16) %z = [zavesch/zflag zavedir/zflag zz Tperiod]; %z = [t; zsch]; % % % Plot of numerically iterated version of check code above vs. known % % solution. % %plot(x,real(ifft(v0)),'--b',x,real(ifft(m0)),'ok',x,real(exp(i*(T))*sin(x)),'r') % % % Plot of $|u_\gamma^\infty|$, $|n_\gamma^\infty|$ % %plot(x,abs(ifft(v0)),'--b',x,abs(ifft(m0)),'ok'); % %Recording the $\dot{H}^1$ norm of $u$ and the $L^2$ norm of $n$. % if flag > 1e-6 % z1 = dx*(sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); % z2 = dx*(sum(abs(ifft(m0)).^2))^(1/2); % z = [z1; z2; 1]; % else % z1 = dx*(sum(abs(ifft(v0)).^2+ abs(ifft(xi.*v0)).^2))^(1/2); % z2 = dx*(sum(abs(ifft(m0)).^2))^(1/2); % z = [z1; z2; 0]; % end