Split-Step Algorithm

This is a fast method for the solution of parabolic PDE’s, relying on the FFT implementation of the Fourier Transform to speed things up. Here it’s applied to an acoustic/seismic problem, following the development of Kuperman/Jackson in “Imaging of Complex Media with Acoustic and Seismic Waves”. Consider first the Helmholtz equation for a point source at r’,z’:

where G is the Green’s function, K the wavenumber (function of the frequency omega and sound speed c). Assuming azimuthal symmetry, G may be expressed as a product of two functions:

and similarly K, now a product of (constant) K0 and index of refraction n. Substituting into the Helmholtz equation gives two PDE’s:

The first PDE has Bessel functions as solutions; taking the assymptotic outgoing Hankel function solution and substituting it into the second, with the narrow angle approximation (second derivative of psi with respect to r much smaller than first derivative wrt r), one finds for the second (parabolic) PDE:

where chi is the fourier transform of psi. Assuming the variation of n is insignificant, in the wave-space domain the PDE and solution are:

Finally, the inverse FT gives the field as a function of depth (Delta r = r-r0, r0 the boundary value):

The following is a little hack in octave, to demonstrate the solution method:

function [Psi,s]=test_ssa(N,Dr,P0,n,K0)

% to test split step algo on parabolic PDE:
%
% \frac{\partial^2 \psi}{\partial r^2}
% +2iK0\frac{\partial \psi}{\partial r}
% +K0^2(n^2-1)\psi = 0
%
% K0*n^2 = (const.) wave number
% N=marching steps in field
% P0=boundary field data
% Dr=field step size

% wjb 0609

%initialize
m=length(P0); Psi=zeros(N,m); Psi(1,:)=P0;

%Nyquist thm, s (wave space) range

s=fftshift(-m/2:m/2-1)./Dr;

for k=2:N

Psi(k,:)=exp((i*K0/2)*(n^2-1)*Dr).*ifft(exp(-(i*Dr*s.^2)./(2*K0)).*fft(Psi(k-1,:)));

end

Depeche Mode

I may just be hanging on to the last tattered threads of my genX youth, but the new album is sublime, particularly tracks 1,4,8,10

Custom Electric Potential

I would like an electric potential with specific qualities, namely, quadratic in two dimensions. On the one hand, I could try and do a little thought experiment, as to which boundary gives the desired properties. Realistically though on the length scale I need it, the odds of producing the boundary are zilch. Much easier in terms of the multipole expansion. With this in mind, and knowing that copper wire along a cylinder makes for an easy setup, and that at least four terms are required for a quadratic potential, the following should work: V(x=0,z=a)=k, V(x=0,z=-a)=k, V(x=-a,z=0)=-k, V(x=a,z=0)=-k, with each wire along y. Assuming long wires, the potential for all four in the x/z plane is:

To confirm this has the desired behavior near the origin (and hopefully over some reasonable distance), I turn to the Taylor multi-variable expansion. Example partial derivatives:

Substituting into the expansion:

does indeed reveal that V ~ (2/{a*a})[z*z-x*x], as desired: