#Wave equation superluminal update
% A wave travelling on a string with % fixed endpoints function main () % KSmrq's colors red = blue = / 256 green = / 256 yellow = / 256 white = 0.99 * % length of the string and the grid L = 5 N = 151 X = linspace ( 0, L, N ) h = X ( 2 ) - X ( 1 ) % space grid size c = 0.5 % speed of the wave tau = 0.25 * h / c % time grid size K = 5 % steepness of the bump S = 0 % shift the wave f = inline ( 'exp(-K*(x-S).^2)', 'x', 'S', 'K' ) % a gaussian as an initial wave df = inline ( '-2*K*(x-S).*exp(-K*(x-S).^2)', 'x', 'S', 'K' ) % derivative of f % wave at time 0 and tau U0 = 0 * f ( X, S, K ) U1 = U0 - 2 * tau * c * df ( X, S, K ) U = 0 * U0 % current U Big = 10000 Ut = zeros ( Big, N ) Ut ( 1, :) = U0 Ut ( 2, :) = U1 % hack to capture the first period of the wave min_k = 2 * N k_old = min_k turn_on = 0 for j = 3 : Big last_j = j % fixed end points U ( 1 )= 0 U ( N )= 0 % finite difference discretization in time for i = 2 :( N - 1 ) U ( i ) = ( c * tau / h )^ 2 * ( U1 ( i + 1 ) - 2 * U1 ( i ) + U1 ( i - 1 )) + 2 * U1 ( i ) - U0 ( i ) end Ut ( j, :) = U % update info, for the next iteration U0 = U1 U1 = U % hack to capture the first period of the wave k = find ( abs ( U ) = max ( abs ( U )) ) k = k ( 1 ) if k > N / 2 turn_on = 1 end min_k = min ( min_k, k_old ) if k > min_k & min_k = k_old & turn_on = 1 break end k_old = k end % truncate to the first period last_j = last_j - 1 Ut = Ut ( 1 : last_j, :) % shift the wave by a certain amount shift = floor ( last_j / 4 ) Vt = Ut Ut (( last_j - shift + 1 ): last_j, :) = Vt ( 1 : shift, :) Ut ( 1 :( last_j - shift ), :) = Vt (( shift + 1 ): last_j, :) num_frames = 100 spacing = floor ( last_j / num_frames ) % plot the wave for j = 1 :( last_j - spacing + 1 ) U = Ut ( j, :) if rem ( j, spacing ) = 1 figure ( 1 ) clf hold on axis equal axis off lw = 3 % linewidth plot ( X, U, 'color', red, 'linewidth', lw ) % plot the ends of the string small_rad = 0.06 ball ( 0, 0, small_rad, red ) ball ( L, 0, small_rad, red ) % size of the window ys = 1.1 axis () % small markers to keep the bounding box fixed when saving to eps plot ( - small_rad, ys, '*', 'color', white ) plot ( L + small_rad, - ys, '*', 'color', white ) frame_no = floor ( j / spacing ) + 1 frame = sprintf ( 'Frame%d.eps', 1000 + frame_no ) disp ( frame ) saveas ( gcf, frame, 'psc2' ) end end function ball (x, y, radius, color ) % draw a ball of given uniform color Theta = 0 : 0.1 : 2 * pi X = radius * cos ( Theta ) + x Y = radius * sin ( Theta ) + y H = fill ( X, Y, color ) set ( H, 'EdgeColor', color ) % The gif image was creating with the command % convert -antialias -loop 10000 -delay 15 -compress LZW Frame10* Movie.gif