% Usage: [u, v] = pset3prob2(x2, T, b, c, lambda, sigma2, stfreq, stcen, sxloc) % % Solves the leap-frog wave equation for periodic boundaries with a source: % exp(-(t-stcen)^2) * sin(stfreq * (t-stcen)) % located at x=sxloc. (For the problem set, stcen = stfreq = 5, and sxloc = 10.) % % x2 is the grid of x points in the unit cell, at *twice* the resolution. % That is, it includes *both* the u and v staggered-grid points. % (For the problem set, x2 = [0:0.05:19.95].) % % T is the time to run for. We start with initial conditions u = v = 0. % b and c are the constants in the u and v equations, and lambda = dt/dx. % % sigma2 is the PML periodicity sigma(x) at the x2 points. % % Returns u and v at time T. u is at x2(1:2:end) and v is at x2(2:2:end). function [u, v] = pset3prob2(x2, T, b, c, lambda, sigma2, stfreq, stcen, sxloc) xu = x2(1:2:end); xv = x2(2:2:end); u = zeros(1,length(xu)); v = zeros(1,length(xv)); % note: we must have length(xu) == length(xv) sigu = sigma2(1:2:end); sigv = sigma2(2:2:end); % h = semilogy(xu, abs(u)); % axis([min(x2) max(x2) 1e-12 1]); h = plot(xu, u); axis([min(x2) max(x2) -1 1]); set(h,'EraseMode', 'xor'); dx = (max(xu) - min(xu)) / (length(xu) - 1); dt = dx * lambda; t = 0; sx = (abs(xu - sxloc) == min(abs(xu - sxloc))) / dx; sx0 = find(sx ~= 0); sx0 = sx0(2:end); sx(sx0) = 0; while (t < T) t = t + dt; stval = exp(-(t-stcen)^2) * sin(stfreq*(t-stcen)); u = (u .* (1 - dx*sigu/2) + lambda * b * (v - [v(end),v(1:end-1)]) + dt * sx * stval) ./ (1 + dx*sigu/2); v = (v .* (1 - dx*sigv/2) + lambda * c * ([u(2:end),u(1)] - u)) ./ (1 + dx*sigv/2); set(h, 'XData', xu, 'YData', u); % set(h, 'XData', xu, 'YData', abs(u)); drawnow; end