Hi Chao,
On point (1): what you think is a perfect delta function, because you only see a single 1 and many zeros, can be thought of as a discrete representation of a continuous signal, i.e. you are only looking at a sampled version of the underlying continuous signal. The underlying continuous signal depends on the basis functions you choose (i.e. how you choose to interpolate between the samples) and for a Fourier basis, the underlying function that fits these samples is a bandlimited, oscillatory function. Run this code and you will see what I mean.
% create the computational grid
Nx = 64; % number of grid points in the x (row) direction
dx = 1; % grid point spacing in the x direction [m]
x = (-Nx/2:Nx/2-1)*dx;
% define a single source pixel (numerical delta function)
delta = zeros(Nx,1);
delta(Nx/2+1) = 1;
% calculate underlying bandlimited delta function using a much finer mesh
dx_fine = 0.001;
x_fine = min(x):dx_fine:max(x);
underlying_func = (1/Nx)*cot(x_fine*pi/(dx*Nx)).*sin(x_fine*pi/dx);
underlying_func(x_fine==0) = 1;
% plots
figure
subplot(3,1,1)
stem(x,delta,'o')
axis([min(x) max(x) -0.22 1])
title('''discrete'' delta function')
subplot(3,1,2)
plot(x_fine,underlying_func,'k')
axis([min(x) max(x) -0.22 1])
title('underlying bandlimited function')
subplot(3,1,3)
stem(x,delta,'o')
hold on
plot(x_fine,underlying_func,'k')
axis([min(x) max(x) -0.22 1])
title('''discrete'' delta is sampled version of bandlimited function')
(2) Yes. You could think of this as a linear system (as long as you stick to the linear wave equation and not the nonlinear version). You can imagine a matrix equation with a vector of p(x,t) stacked up on the left, a matrix made up by discretising the wave operator, and an initial condition (or a source, depending on how you have written the operator matrix) on the right. You haven't avoided solving the wave equation: this just is a method for solving the wave equation (and probably not a very efficient one). You need to choose some way of discretizing the wave operator to form the system matrix, and that might be using finite differences, finite elements, spectral approaches, or whatever you like!
(3) These are equivalent in the sense that p_0 = \Gamma H is true for photoacoustics (when the optical pulse is short enough), and dp/dt = 0 when u = 0 (which can be seen by combining the mass conservation equation and pressure-density relation).
Hope that helps,
Ben
