function aufgabe9a(N) if (nargin < 1) N = 256; end V = @(x) -cos(x/2); %s = 0.1; psi0 = @(x) exp(-10*(x-0.5).^2) * (20/pi)^(1/4); mu = 1/0.1^2; k = (-N/2):(N/2-1); x = k * 2 * pi / N; w = ones(size(x)) * (x(2) - x(1)); w(1) = w(1)/2; w(end) = w(end)/2; plot(x, V(x), '-k', x, psi0(x), '-r'); xlim([-pi pi]); drawnow; fprintf('norm is: %g\n', w * psi0(x').^2); DK = 1/(2 * mu) * k.^2; tic F = 1/N * exp(-2i*pi/N * k' * k); FI = exp(2i*pi/N * k' * k); toc %tic; %V1 = F * diag(V(x)) * FI; %toc x1 = linspace(-pi, pi, N+1); w1 = ones(size(x1)) * (x1(2) - x1(1)); w1(1) = w1(1)/2; w1(end) = w1(end)/2; tic V2 = 1/(2*pi) * ones(N,1) * (w1 .* V(x1)) .* ... exp(-1i*k'*x1) * transpose(exp(1i*k'*x1)); toc c0 = F * psi0(x'); %c1 = 1/(2*pi) * exp() (w1 * psi0(x))' %f = @(t, y) -1i * (diag(DK) + V1) * y; %tic %[t1, ct1] = ode45(f, [0 100], c0); %toc %tic %[t2, ct2] = ode15s(f, [0 100], c0); %toc k = [0:(N/2-1) (-N/2):-1]; DK = 1/(2 * mu) * k.^2; c0 = fft(psi0(x')); f = @(t, y) -1i * (DK' .* y + fft(V(x1(2:end))' .* ifft(y))); fprintf('V=cos(x), ode45 '); tic [t3, ct3] = ode45(f, [0 100], c0); toc view2D(1, ct3, t3, x); title('ode45'); fprintf('V=cos(x), ode15s '); tic [t4, ct4] = ode15s(f, [0 100], c0); toc view2D(2, ct4, t4, x); title('ode15s'); h = @(x, a) 0.5*(sign(x-a)+1); hm1 = @(x, a) -0.5*(sign(x-a)-1); %Vs = @(x) hm1(x,-5).*((x+5).^2 + 2) + h(x,-5).*h(x,0).*ceil(-2/5*x) ... % + (max(x,5)-5).^2; %V = @(x) Vs(x*pi/10); V = @(x) abs(x); %plot(x, V(x)); c0 = fft(psi0(x')); f = @(t, y) -1i * (DK' .* y + fft(V(x1(2:end))' .* ifft(y))); fprintf('V=|x|, ode45s '); tic [t5, ct5] = ode45(f, [0 100], c0); toc view2D(3, ct5, t5, x); title('V(x) = |x|'); figure(4); clf; Z = ifft(transpose(ct5)); NX = length(x); IX = 1:ceil(NX/20):NX; NT = length(t5); IT = 1:ceil(NT/50):NT; [X T] = ndgrid(x(IX), t5(IT)); quiver(X, T, real(Z(IX,IT)), imag(Z(IX,IT))); return function view2D(fig, ct, t, x) [X T] = ndgrid(x, t); Z = ifft(transpose(ct)); C = thaller(Z, 0.5); hf = figure(fig); clf; contour3(X, T, abs(Z), 4, '-w'); % Konturlinien hs = surface(X, T, abs(Z), hsv2rgb(C), 'EdgeColor', 'none', ... 'FaceColor', 'interp'); % Oberflaeche view(-20, 60); % netter Blickwinkel return function C = thaller(Z, R) [Nx, Ny] = size(Z); % Farbtabelle gemaess Thaller h = mod(-angle(Z) / (2*pi), 1); t = pi - 2*atan2(abs(Z), R); l = 1 - t / pi; % HSV Tupel C = ones(Nx, Ny, 3); C(:,:,1) = h; C(:,:,3) = l; return