Summary week 7

1
Summary week 7

• Fast Fourier transform
• Example with gravitational waves, and exercise
• Huge enhancement signal/noise for periodic
  signals
• Wavelets
• Linear operation like FFT
• Daubechy set
• Smoothing operator running weighted average
• Cancels constant/linear terms
• Difference operator
• Enhances the local structures
• Wavelets are located in both domains
• Wavelets have discontinuous derivatives
• Reduction of information

2
wavelet transform

• What do these wavelets look like?
• perform DAUB4 and DAUB20 on unit vectors

3
wavelets

• cusps non-smooth derivatives
• use of wavelets
• 1) to get strength of located events (e.g. a
  signal in a measurement stream)
• 2) to render matrices sparse.

4
Traveling salesman problem

• Energy function total distance of path
• Steps
• Moving 1 city (6 distances change)

5
Traveling salesman problem

• Energy function total distance of path
• Steps
• Moving 1 city (6 distances change)
• in principle all paths can be built from this
  step
• Swapping 2 cities (8 distance changes)

6
Traveling salesman problem

• Moving a stretch of k1 cities
• 6 distances to calculate total path lenght, 2
  additional distances if reversal of order is
  included

7
Traveling salesman

• Make choice of design
• I chose Class to read in cities
• IO, distance between 2 cities, index
• Distance I went to polar coordinates and use
• Structures to keep track of different paths
  between cities
• I implemented Vector of integers, indices of
  cities
• Array containing distances between cities i and j
• Method to calculate distance for a given path
• Methods to calculate Delta-E for a proposed
  pathchange (move, swap, move/reverse stretch)
• Store current distance, calculate change

8
Traveling salesman

• Make choice of design
• Methods to execute path changes
• Methods for the random-number generation
• I used NRran1 and NRexpdev, that is expensive
• (methods to keep track of statistics)
• In main
• Control value of Delta-E
• Control which methods are used with what
  parameters (lenght of the stretch of cities,
  reversal, start configurations, number of trials)
• Keep track of shortest path encountered
• Verify end result (millions of changes
  recalculation of total path is the same within
  micrometers).
• Next analyze system path lenghts, parameters

9
Traveling Salesman

• Start configuration
• Random average distance about 1.2-1.3 Gm
• Short take random city, take closest neighbour
• Average distance 165 Mm

How to calculate amount Of possible
configurations?
Overflows on 64-bit machine Take logarithms of i,
add them, take exponent of remainder
10
Traveling Salesman

• Explore strategies values of sigma, number of
  trials, methods, path lenghts
• I started with 10 start configurations
• Ntrial 20000, sigma 5000 km
• If after 20000 tries no shorter path is found,
• Reduce sigma with factor 3
• Increase amount of trials with factor 2
• Until sigma lt 1 km
• Around 1 million different configurations, takes
  a fraction of a second
• How many configurations to study?
• Keep track of improvement as a function of amount
  of trials

11
Parameter determination
method. sigma mean av move1 37 150.6 6.37 mov
e1 1-12 150.2 6.18 swap 37 202.7 8.92 swap 1-12
202.3 8.754 stretch 37 155.0 4.23 stretch 12 15
1.3 4.39 stretch 4 146.4 2.48 stretch 1 143.6 1.
93 all 37 134.1 1.672 all 12 133.6 1.592 all 4
133.0 1.571 all 1 132.5 1.803
All methods together shortest paths Gain from
sigma 12-1 1000 km Computation time loss of
factor 3
12
Search strategy

• Shortest path is found when all 3 methods are
  used together
• Sigma 5000 too large (after 20000 trials paths
  still around 700,000 km)
• Start with sigma 1000, path lenghts up to 60
• For sigmas smaller than 12 not much improvement
• Trials went up to 640000 failures
• 3 times a shorter path found, average lt1000 km
  shorter
• Path lenghts reductions for swapping up to 60
  cities in initial stages (sigma gt100)
• Reductions up to 8 cities sometimes at small
  sigma (10 km)
• My chosen strategy
• start with sigma of 1000, divide by 2.5 each
  refinement
• Start with 20000 trials, multiply by 1.5 each
  refinement
• Stop at sigma lt1, if distance smaller than
  129.500 km, try more trials at sigma of 10, 4,
  1.5 (triples computation time for shortest 5 of
  paths)

13
Final results

• 1500 start configurations in 12 minutes
• 20 short start paths pick closest city
• Shortest distance 127096 km
• (better strategy lt126000, worse strategy
  gt129000)
• Reduction of pathlenghts

method short start path sigma mean RMS 1.6 132.2
1819 4 132.2 1802 10 132.5 1827 25 133.2 1937 6
4 134.8 1876 160 145.8 4098 400 164.3 10580 1000 1
64.3 10580
method all sigma mean RMS 1.6 131.8 1693 4 131.9
1714 10 132.2 1776 25 132.8 1866 64 134.4 2147 160
143.7 2827 1000 310.4 8014
14
Final Results
White sigma1 Red sigma4 Blue
sigma10 Turquoise sigma25 Yellow
sigma64 Pink sigma160
15
Final results

• 3 distances lt 128000 km in 800 seconds
• 79 configurations lt 129500 km

16
Integration of functions (Quadrature)

• intuitive expectation integral is much harder to
  calculate than differential
• basic task calculate
• chapter 16 integration of differential
  equations.
• if possible, use this differential form when the
  function is concentrated in sharp peaks or when
  the shape is strongly scale-dependent.
• Already discussed Chebyshev integration
• Fourier integration also methods available

17
Classical formulas
18
Trapezoidal Rule

• zeroth- order function replaced by sum

19
Trapezoidal rule, Simpsons rule

• The trapezoidal rule is correct for first-order
  polynomials. It considers 2 points, ½(f0f1),
  so it is called a two-point function
• With more than 2 points, one can calculate an
  integral between these points to higher order.
  The three-point formula is correct to order 2
  (and due to a cancellation, also to order 3)
  Simpsons rule
• The four-point formula is also correct up to
  order three.

20
Classical formulas

• The five-point formula is correct to fifth order
  Bodes formula
• factors derive from writing down the equations
  with arbitrary factors p,q,r,s,... calculate the
  integral for functions 1, x, x2, x3, ... and
  solving the coupled equations.

21
Extended formulas (Closed)

• calculate N times for intervals 0-1, 1-2,...
• Example 5 points calculation of sin(x) between
  0,pi and ex and sqrt(x) between 0,2, with
  trapezoidal, Simpsons, and Bodes rule

22
Extended formulas (Closed)

• calculate N times for intervals 0-1, 1-2,...
• Example 5 points calculation of sin(x) between
  0,pi and ex and sqrt(x) between 0,2, with
  trapezoidal, Simpsons, and Bodes rule

23
Extended formulas (Closed)

• calculate N times for intervals 0-1, 1-2,...
• Example 5 points calculation of sin(x) between
  0,pi and ex and sqrt(x) between 0,2, with
  trapezoidal, Simpsons, and Bodes rule

24
Extended formulas

• Precision is gained, just by weighting the
  function values differently!
• 1 order of precision obtained, just from relative
  weight of end point ( both for Simpson and for
  trapezoidal rule).
• Open formulas

25
Extended formulas (open)

• construct them by using the closed formula for
  the interior part and adding the extrapolated
  open step for the ends.
• end steps done only once. Use 1 order lower
  extrapolation than for interior steps.
• open formulas 1 order lower precision overall.

26
Elementary algorithms

• Starting point extended trapezoidal rule
• one can add new points in between without loosing
  the previous work.
• error extrapolation is entirely even in powers of
  1/N

27
Euler-MacLaurin Summation formula

• Only even terms enter in this series.
• calculate integral with N steps (SN), then with
  2N steps (S2N)

28
Romberg integration

• The leading error in the second step will be ¼ of
  the previous step. It can be cancelled by
• The next order error is fourth order in 1/N. This
  procedure, applied once, is exactly equivalent to
  Simpsons rule.
• Romberg integration use k successive refinements
  to eliminate errors up to O(1/N2k).
• Much more precise than trapezoid rule or Simpsons
  rule.

29
Improper integrals

• Improper
• finite value, but cannot be calculated at
  boundary (e.g.) (sinx/x at 0)
• one of the limits is infinity
• integrable singularity (e.g. 1/sqrt(x) at 0)
• at boundary
• at known place in interval a,b
• at unknown place
• How to solve we need an algorithm like the
  Romberg integration, but for open formula.
  Extended midpoint rule also only has even errors

30
Open integration

• It is not possible to double the number of points
  by adding intermediate points, but you can triple
  them.
• the new sum is S(9S3N-SN)/8.
• routine qromo does Romberg-like integration for
  open intervals.
• infinite boundary-gt change of variables
• if the range a-b contains 0, split the integral
  in two parts. (a-1, 1-inf)

31
Gaussian quadrature

• Freedom to choose coefficients (Romberg
  integration) leads to much faster convergence.
• Freedom to choose abscissas twice as much
  freedom
• you can arrive at higher-order convergence
  faster, when integrand is smooth
• you can make the integral exact for a class of
  functions of polynomialsweight function W

32
Gauss-Chebyshev

• e.g. integral

33
Gaussian quadrature

• Chebyshev W(x) 1/sqrt(1-x2)
• W1 tabulated weights qgaus (10 points)
• How to calculate weights and abscissa?
• orthonormal set ltfggt 0 if f not equal g,
  ltffgt1
• orthonormal sets can be found by recurrence

34
Gaussian quadrature

• The polynomials pj have exactly j distinct roots.
  The roots of pj interleave the roots of pj1.
• The abscissa in Gaussian quadrature are formed by
  the roots (Like in Chebyshev)
• The weights can be calculated by e.g.

35
Gaussian quadrature

• examples
• routines are provided in the book.

36
Multi-dimensional integrals

• number of function evaluations grow like the
  power of the dimension (e.g. typically 30 for
  1-dim, 30,000 for 3 dim)
• boundary may be complicated
• Try always to reduce it to lower dimension.
• (e.g. spherical symmetry polar coordinates)
• complicated boundary resort to Monte Carlo
  techniques.
• care when function is strongly peaked
• if boundary is simple, Gaussian quadratures may
  give a relatively fast answer.

37
Multidimensional integration