previous up next print clean
Next: Side-boundary analysis Up: WAVEMOVIE PROGRAM Previous: Time-domain analysis

Internals of the film-loop program

The differential equation solved by the program is  
 \begin{displaymath}
{\partial P \over \partial z } \eq 
{ i \, \omega \over v(x,...
 ...er - \, i \, \omega \, 2 } \ 
{\partial^2 P \over \partial x^2}\end{displaymath} (53)
For each $\Delta\,z$-step the calculation is done in two stages. The first stage is to solve  
 \begin{displaymath}
{\partial P \over \partial z} \ \eq 
{v \over - \, i \, \omega \, 2} \ {\partial^2 P \over \partial x^2}\end{displaymath} (54)
Using the Crank-Nicolson differencing method this becomes
\begin{displaymath}
{ p_{{z+1}}^x - p_z^x \over \Delta z }
\eq
{v \over -i\omega...
 ...p_{{z+1}}^x + 
p_{{z+1}}^{{x-1}} \over 2\, \Delta x^2 } \right)\end{displaymath} (55)
Absorb all the constants into one and define  
 \begin{displaymath}
\alpha \ \eq \ {v \, \Delta z \over -\,i\,\omega\,4\, \Delta x^2 }\end{displaymath} (56)
getting
\begin{displaymath}
p_{{z+1}}^x \,-\, p_z^x
\eq
\alpha \, \left[
\, ( \, p_z^{{x...
 ...+1}}\,-\, 2 \, p_{{z+1}}^x \,+\, p_{{z+1}}^{{x-1}} )
\, \right]\end{displaymath} (57)
Bring the unknowns to the left:  
 \begin{displaymath}
- \alpha p_{{z+1}}^{{x+1}} \,+\,
( 1+2 \alpha ) p_{{z+1}}^x ...
 ...p_z^{{x+1}} \,+\,
(1-2 \alpha ) p_z^x \,+\,
 \alpha p_z^{{x-1}}\end{displaymath} (58)
We will solve this as we solved equations (37) and (38). The second stage is to solve the equation  
 \begin{displaymath}
{\partial P \over \partial z} \ \eq \ {i\,\omega \over v }\ \ P\end{displaymath} (59)
analytically by  
 \begin{displaymath}
P(z+ \Delta z) \ \eq \ P(z) \ e^{i \Delta z \ \omega / v }\end{displaymath} (60)

By alternating between (58) and (60), which are derived from (54) and (59), the program solves (53) by a splitting method. Formal justification of the splitting method follows in chapter [*]. The program uses the tridiagonal solver discussed earlier, except the version here, ctris(), has all the real variables changed to complex.

Figure 5 shows a change of initial conditions where the incoming wave on the top frame is defined to be an impulse, namely, p(x,z=0) = $(\cdots,0,0,1,0,0,\cdots)$.The result is alarmingly noisy. What is happening is that for any frequencies anywhere near the Nyquist frequency, the difference equation departs from the differential equation that it should mimic. This problem is addressed, analyzed, and ameliorated in chapter [*]. For now, the best thing to do is to avoid sharp corners in the initial wave field.

 
Mcompart
Figure 5
Observe and describe various computational artifacts by testing the program using a point source at (x,z) = (xmax/2,0). Such a source is rich in the high spatial frequencies for which difference equations do not mimic their differential counterparts.

Mcompart
view burn build edit restore


previous up next print clean
Next: Side-boundary analysis Up: WAVEMOVIE PROGRAM Previous: Time-domain analysis
Stanford Exploration Project
10/31/1997