Processing math: 0%

Approximate Density

\def \half {\frac{1}{2}} \def \hieruit {\quad \Longrightarrow \quad} \def \slechts {\quad \Longleftrightarrow \quad} \def \OF {\quad \mbox{or} \quad} \def \MET {\quad \mbox{with} \quad} The reverse question. Suppose that the sense P(x) is given, how can we find then the associated density of points \{x_k\} ? A crude approximation is to take the value P(x_k) at the position x_k and construct the next point by demanding that the rectangle with height P(x_k) and width (x_{k+1}-x_k) has an area equal to one: P(x_k) (x_{k+1}-x_k) = 1 \hieruit x_{k+1} = x_k + \frac{1}{P(x_k)} Herewith it is assumed that P(x) is positive everywhere. But, as we have seen with the guitar example in the previous section, it is also possible that P(x) is negative everywhere, giving rise to a monotone decreasing sequence. In that case, the x-axis shall be traversed in the opposite direction, giving: P(x_k) (x_k-x_{k-1}) = 1 \hieruit x_{k-1} = x_k - \frac{1}{P(x_k)} Anyway, a much better job can be done as follows. Consider a linear approximation of the sense in the neighbourhood of x_k: P_{k+1} = P(x_k) + P'(x_k) (x_{k+1} - x_k) Where it is noticed that the estimate P_{k+1} may be almost, but not exactly equal to P(x_{k+1}). The density may be approximated with help of the trapezium rule as well: \half \left[ P_{k+1} + P(x_k) \right] \left[x_{k+1} - x_k\right] = 1 Substitution of the first formula into the second gives: \half \left[ P(x_k) + P'(x_k) (x_{k+1} - x_k) + P(x_k) \right] \left[x_{k+1} - x_k\right] = 1 \hieruit \left[ \half P'(x_k) \right] (x_{k+1} - x_k)^2 + P(x_k)(x_{k+1} - x_k) - 1 = 0 A quadratic equation. Solve for (x_{k+1} - x_k) : x_{k+1} - x_k = \frac{ - P(x_k) \pm \sqrt{P^2(x_k) + 2 P'(x_k)}}{P'(x_k)} Take as an example the linear density: P(x) = \frac{1}{\Delta^2}\left(x+\frac{\Delta}{2}\right) \hieruit P'(x)= \frac{1}{\Delta^2} \hieruit \\ x_{k+1} - x_k = \frac{ - (x_k+\Delta/2)/\Delta^2 \pm \sqrt{\left[(x_k+\Delta/2)/\Delta^2\right]^2+ 2/\Delta^2}}{1/\Delta^2} \hieruit \\ x_{k+1} = \pm \sqrt{(x_k+\Delta/2)^2 + 2\Delta^2}-\Delta/2 Compare with a previously found result to establish that the sign of \pm must be plus: x_k = \left( \sqrt{1/4 + 2 k}-1/2 \right) \Delta \hieruit (x_k+\Delta/2)^2 = \Delta^2(1/4 + 2 k)^2 \hieruit \\ x_{k+1} = \pm \Delta\sqrt{(1/4 + 2 k) + 2}-\Delta/2 = \left( \sqrt{1/4 + 2 (k+1)}-1/2 \right) \Delta It should be no surprise that the outcome for a linear density is not approximate but exact.
Anyway, quite in general, we still have: x_{k+1} - x_k = \frac{\sqrt{P^2(x_k) + 2 P'(x_k)} - P(x_k)}{P'(x_k)} The numerator is of the form (a-b). Now multiply with (a+b)/(a+b), leading to values of (a^2-b^2)/(a+b) which are likely more stable numerically: x_{k+1} - x_k = \frac{2 P'(x_k)/P'(x_k)}{\sqrt{P^2(x_k) + 2 P'(x_k)} + P(x_k)} \\ \hieruit x_{k+1} = x_k + \frac{1}{\sqrt{\left[P(x_k)/2\right]^2 + P'(x_k)/2} + P(x_k)/2} The plus sign (out of \pm\sqrt{}) also comes from the fact that, for P'(x) = 0, the next point has to be found by the (no longer) crude approximation: x_{k+1} = x_k + 1/P(x_k). It is noticed, in addition, that the expression under the square root should not become negative. Meaning that, if P'(x_k) < 0: | P'(x_k) | < \half \left[ P(x_k) \right]^2 It can be argued, though, that the above condition must also hold for positive values of P'(x_k). The reason is that the sequence x_k could be traversed in the opposite direction as well: \left. \begin{array}{l} f_{k-1} = P(x_k) + P'(x_k) (x_{k-1} - x_k) \\ \half \left[ P(x_k) + f_{k-1} \right] \left[x_{k-1} - x_k\right] = -1 \end{array} \right\} \hieruit \left[ \half P'(x_k) \right] (x_{k-1} - x_k)^2 + P(x_k)(x_{k-1} - x_k) + 1 = 0 \left[ \half P'(x_k) \right] (x_{k+1} - x_k)^2 + P(x_k)(x_{k+1} - x_k) - 1 = 0 Thus (k+1) changes into (k-1) means that +1 changes into -1. And that's the only difference between traversing in a positive or in a negative direction. A bit of further thinking reveals then that the negative end-result must be: x_{k-1} = x_k + \frac{1}{\sqrt{\left[P(x_k)/2\right]^2 - P'(x_k)/2} + P(x_k)/2} And therefore, again, the abovementioned condition must hold: | P'(x_k) | < \half \left[ P(x_k) \right]^2 Take for example the linear density, where P(x) = (x/\Delta + \half)/\Delta, P'(x) = 1/\Delta^2. Then our condition translates into: | P'(x_k) | < \half \left[ P(x_k) \right]^2 \slechts 1/\Delta^2 < \half \left[ (x_k/\Delta + \half)/\Delta \right]^2 \slechts \left[ (x_k/\Delta + \half) \right]^2 > 2 \slechts x_k/\Delta > - \half + \sqrt{2} \OF x_k/\Delta < - \half - \sqrt{2} Only the non-negative solution must be kept. We conclude that the algoritm for constructing the set \{x_k\} from the linear sense P(x) = (x/\Delta + \half)/\Delta will certainly work for x_k > (\sqrt{2} - 1/2) \Delta \approx 0.9142135 \Delta , let's say x_k > \Delta, which is certainly the case for k > 0.
Take as an example the sense that is an approximation for the longitudinal wave with small amplitude, eventually without the latter restriction: P(x_k) = \frac{1}{\Delta}\left[1 - A\frac{2\pi}{\lambda}\cos\left(\frac{2\pi}{\lambda} x_k\right)\right] \hieruit \\ P'(x_k) = \frac{A}{\Delta}\left(\frac{2\pi}{\lambda}\right)^2\sin\left(\frac{2\pi}{\lambda} x_k\right) From which the iteration scheme follows. The result is:

Accompanying software for making the pictures is disclosed as well.