tnotes.rst - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tnotes.rst (13063B)
---
1 # vim: set tw=80:noautoindent
2
3 For Jacobi exercise: Use atan2().
4
5 11-4: Introduction
6 #####
7 Deadline for exercises: ~2 weeks after assignment
8
9 All exercises must be done on Lifa:
10 lifa.phys.au.dk
11 Ethernet: vlifa01 or vlifa02
12 You need RSA keys to login from outside net.
13 Login from inside: NFIT credentials.
14
15 Alternative server: molveno
16 Dedicated for numerical course.
17
18 Lifa will have approx. 5 year old software, molveno is updated.
19
20 All phys.au.dk servers have NFS shared home folders.
21
22 Dmitri answers:
23 http://www.phys.au.dk/~fedorov/numeric/
24
25 Bash: Last command starting with e.g. 'v': !v
26
27 Exercises are weighted 75% and the exam 25%. You do need to complete at least
28 51% of the exam.
29
30 02: >51%
31 4: >60%
32 7: >80%
33 10: >90%
34 12: 100%
35
36 16-4: Interpolation
37 ####
38 Makefiles_:
39 For all project, use makefiles, the project is evaluated by `make clean; make`.
40
41 General structure of Makefile:
42 Definitions of variables
43 all: my_target
44 Target(s): Dependenc(y,ies)
45 <-tab-> Command(s)
46
47 Strings are denoted without decoration, variables as $(VAR).
48 To use the CC system linker to link C++ objects, add libraries:
49 LDLIBS += -lstdc++
50 If you use e.g. g++ as the linker, the above command is not needed.
51
52 If an object file is required as a dependency, the Makefile will issue a CC
53 command, compiling the source code file with the same basename. `make -p | less`
54 will show all build rules, even build in.
55
56 Interpolation_:
57 When you have a discrete set of datapoints, and you want a function that
58 describes the points.
59
60 If a function is analytical and continuous, and an infinitely large table of
61 datapoints in an interval, the complete analytical function can be found from a
62 taylor-series of _all_ derivatives in the interval.
63
64 k-Polynomial: k+1 unknowns (variables). High polynomials tend to oscillate.
65
66 Cubic interpolation: The interval between points can be described by a function
67 of three polynomials.
68
69 Spline interpolation: Also made of polynomials. First order spline (a0+x*a1),
70 simply connects the data points linearly. The splines for each inter-datapoint
71 interval must have the same values at the data points. The derivates are
72 discontinuous.
73
74 Solving linear splines:
75 Datapoints: {x_i,y_i}, i = 1,...,n
76 Splines (n-1): S_i(x) = a_i + b_i (x - x_i)
77 S_i(x_i) = y_i (n equations)
78 S_i(x_{i+1}) = S_{i+1}(x_i) (n-2 equations (inner points))
79 Unkowns: a,b: 2n-2 variables. Solution:
80 => a_i = y_i
81 S_i(x) = y_i + b_i (x-x_i)
82 S_i(x_{i+1}) = y_i + b_i (x_{i+1} - x_i)
83 => b_i = (y_{i+1} - y_i) / (x_{i+1} - x_i)
84
85 Start of by searching which interval x is located in. Optionally, remember the
86 interval, as the next value will probably lie in the same interval.
87
88 Binary search: Is x between x_1 and x_n? If not, error: Extrapolation.
89 Continuously divide the interval into two, and find out if the x is in the
90 larger or smaller part.
91 DO A BINARY SEARCH.
92
93 Second order spline:
94 S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2
95 Unknowns: 3(n-1) = 3*n-3
96 Equations: 3*n-4
97 The derivatives are also continuous.
98
99 Solution:
100 a_i = y_i
101 \delta x = (x-x_i)
102 y_i + b_i \delta x_i + c_i \delta x_i^2 = y_{i+1}
103 b_i + 2 c_i \delta x_i = b_{i+1}
104 Suppose you know b_i, you can find c_i. From that you can find b_{i+1}, and in
105 turn, c_{i+1}. Through recursion you find all unknowns by stepping forward.
106 The backwards solution can be found from the last data-point (y_n) by solving
107 the two equations with c_{n-1} and b_{n-1} as the two unkowns.
108
109 Symmetry can be used as an extra condition.
110
111
112 18-4:
113 ##############
114 Interpolation exercise, plotting optionally in gnuplot, or graph (from
115 plotutils):
116 graph -T png points.dat > plot.png
117 In makefile, example:
118 plot.png: points.dat
119 graph --display-type png $^ > $@
120 Each dataset in points.dat needs a header, e.g. # m=1, S=0
121
122 lspline.cc: Linear spline
123 qspline.cc: Quadratic spline
124 cspline.cc: Cubic spline
125
126 Linear spline:
127 S(x) = S_i(x) it x in [x_i,x_(i+1)]
128 S_i(x) = a_i + b_i (x-x_i)
129 S_i(x) = y_i + (\delta y_i)/(\delta x_i) (x-x_i)
130 S_i(x_i) = y_i
131 \delta y_i = y_(i+1) - y_i
132 S_i (x_(i+1)) = y_(i+1)
133
134 In C++:
135 std::vector<double> x,y,p;
136 Maybe typecasting? Could be fun.
137
138 Procedural programming:
139 --
140 struct lspline {
141 int n;
142 vector<double> x,y,p;
143 };
144
145 Make a function:
146 struct lspline* construct_lspline(vector<double>&x, vector<double>&y)
147 double evaluate_lspline(struct lspline * asdf, double z)
148 free_lspline();
149
150 Object-oriented programming:
151 --
152 If you want to take the structural approach, you keep the functions and
153 structures seperate. If you take a OO approach, put the functions inside the
154 structure (or class):
155 struct lspline {
156 int n;
157 vector<double> x,y,p;
158 lspline(...,..); // Constructor, same name as struct
159 double eval(double z);
160 };
161 struct lspline ms(x,y);
162 ms.eval(z);
163
164 See Dmitri's cubic spline example which uses OO.
165
166 Functional programming (in C++11), compile with -std=c++0x:
167 --
168 The functions can return functions:
169
170 #include <functional>
171 using namespace std;
172
173 function<double(double)> flspline (vector<double>&x, vector<double>&y);
174
175 auto my_spline = flspline(x,y);
176 my_spline(5,0);
177
178
179 System of linear equations:
180 -------
181 A*x = b
182
183 A_i,j x_j = b_i
184
185 Solve by finding A^(-1): x = A^(-1) * b
186 Numerically, you calculate the inverse by solving Ax=b.
187 We will assume that the matrixes are not singular.
188
189 Turn the system into a triangular form.
190 The main diagonal is non-zero, all lower values are 0, and upper values are
191 denoted T_nn.
192 T_nn * x_n = b_n => x_n = 1/T_nn * b_nn
193 T_(n-1,n-1) x_(n-1) + T_(n-1,n) x_n = b_(n-1)
194 T_ii x_i + sum^n_(j=i+1) T_(i,j) x_j = b_i
195 x_i = 1/T_ii (b_i - sum^n_(j=i+1) T_ij, x_j)
196
197 The simplest triangulation is by Gauss elimination. Numerically, the simplest
198 method is LU decomposition (Lower Upper).
199 A = LU, where L=lower triangle, U=upper triangle.
200 n^2 equations.
201 L and U contain "(n^2 - n)/2 + n" elements.
202 L+U = n^2/2 + n/2 = (n(n+1))/2
203
204 The diagonals in L are all equal to 1: L_ii = 1.
205 See Dolittle algorithm in the lecture notes, which with the LU system is the
206 most used, and fastest method for solving a linear equation.
207
208 Ax = b
209 LUx = b, Ux=y
210
211 Another method: QR decomposition: R=Right triangle (equal to U).
212 A = QR
213 Q^T Q = 1 (orthogonal)
214
215 Ax = b
216 QRx = b
217 Rx = Q^T b
218
219 Gram-Schmidt (spelling?) orthogonalization:
220 Consider the columns of your matrix A. Normalize them. Orthogonalize all other
221 columns to the first column.
222
223 Normalizing the column: ||a_1|| = sqrt(a_1 dot a_i)
224 Orthoginalize columns: a_2/||a_1|| -> q_1
225
226 Numerically:
227 for j=2 to m:
228 a_j - dot(dot(q_1,a_j),q_1) -> a_j
229 a_2 -> a_2/||a_2|| = q_2
230
231 Find inverse matrix:
232 A A^-1 = diag(1)
233
234
235
236 30-4: Diagonalization
237 #########################
238
239 Runtime comparison: Do a number of comparisons with different matrix sizes
240 etc.numeric-2012
241
242 Easiest way to diagonalize a matrix: Orthogonal transformation
243 A -> Q^T A Q = D
244 Q matrix can be built with e.g. QR decomposition. Rotation: Computers to cyclic
245 sweep, which is faster than the classic rotation.
246 Cyclic: Zero all elements above the diagonal, and do your rotations until the
247 matrix is diagonal. The matrix is converged if none of the eigenvalues have
248 changed more than machine precision. You will destroy the upper half plus the
249 diagonal. If you store the diagonal in another vector, you can preserve the
250 matrix values.
251 In the beginning, you look for the largest element in each row, and create an
252 index which is a column that keeps the numbers of the largest elements.
253 With an index, you can perform the operation fast. You have to update the index
254 after each operation.
255
256 For very large matrixes, > system memory:
257 Krylov subspace methods: You use only a subspace of the whole space, and
258 diagonalize the matrix in the subspace.
259
260
261 02-5: Ordinary Differential Equations (ODEs)
262 ############################################
263
264 dy/dx = f(x,y)
265
266 Usually coupled equations:
267 dy_1/dx = f_1(x,y_1,y_2)
268 dy_2/dx = f_1(x,y_1,y_2)
269 [y_1, ..., y_n] = y
270 [f_1(x,y), ..., f_n(x,y) = f(x,y)
271 y' = f(x,y)
272 x is usually time. If f has the above form, it is a ODE.
273 Sine function:
274 y'' = -y
275 In this form, it is not a ODE, because it has a second derivative.
276 You can redifine it to be an ODE:
277 => y = u_1 -> u'_1 = u_2
278 y' = u_2 -> u'_2 = -u_1
279
280 Solving ODEs with a starting condition:
281 y' = f(x,y)
282 y(x_0) = y_0
283 a->b (Shape of the function in the interval)
284 You can calculate the derivative at (x_0,y_0). You then make a step towards
285 x_0+h, and apply Euler's method:
286 y' \approx (\delta y)/(\delta x) = (y(x_0+h) - y_0)/h
287 y(x_0+h) \approx y(x_0) + h*f(x_0,y_0)
288 y(x) = y_0 + y'_0(x-x0) + 1/(2!) y''_0(x_0)(x-x_0)^2 + ... (Taylor exp.)
289 You can find the higher-order terms by sampling points around (x_0,y_0).
290 Estimate the new derivative half a step towards h, which is used for the first
291 point. You sample your function, and fit a polynomial. Higher order polynomiums
292 tend to oscillate.
293 When solving ODE's the many sampled points create a tabulated function. If you
294 want to do something further with the function, you interpolate the points by
295 e.g. splines.
296 Only three steps are needed if the equation is a parabola. Other functions are
297 called multi-step methods. You do not want to go to high in the
298 Taylor-expansion, as there lies a danger with the higher order terms
299 oscillating.
300 If the function is changing fast inside an interval, the step size h needs to be
301 small. This is done by a driver.
302 The Runge-kutta is single step, and the most widely used.
303
304 Absolute accuracy (delta) vs. relative accuracy (epsilon) can behave
305 significantly differently when the machine precision is reached.
306 The user usually specifies both, and the accuracy is satisfied, when one of the
307 conditions are met (the tolerance, tau).
308 tau = delta + epsilon*||y||
309
310 Stepper: Must return y(x+h), e (error estimate).
311 Driver: x->x+h. If e is smaller than tau, it accepts the step.
312 The driver finds the size of the next h:
313 h_next = h(tau/e)^power * safety
314 power = 0.25
315 safety = 0.95
316 The driver thus increases the step size if the error was low relative to the
317 tolerance, and vice-versa.
318
319 Runge principle for determining the error:
320 e = C*h^(p+1)
321 You first do a full h step, then two h/2 steps.
322 2*c*(h/2)^(p+1) = c * (h^(p+1))/(2^p)
323 For the error, you calculate y_1 from full step, and then y_1 from two half steps:
324 e = ||y_1(full step) - y_1(2 half steps)||
325 You can also instead use RK1 and RK2, and evaluate the difference between the
326 two for the same step.
327
328 The effectiveness of the driver and stepper is determined by how many times you
329 call the right-hand side of the ODEs.
330
331 Save the x and y values in a C++ vector, and add dynamically using the
332 push_back() function.
333
334
335 07-5: Nonlinear equations and optimization
336 #####
337 Pipe stdin to program:
338 cat input.txt | ./my_prog
339 Redirect stdout to file:
340 ./my_prog > out.txt
341 Redirect stderr to file:
342 ./my_prog 2> err.txt
343
344 Jacobian matrix: Filled with partial derivatives. Used for linearizing the
345 problem. f(x+\delta x) \approx f + J \delta x
346
347 Machine precision: double: 2e-16. It is the largest number where 1 + eps = 1.
348
349 Quasi-Newton methods: Might be exam exercise.
350
351
352 14-5: Numerical integration
353 #####
354 Examination problem: FFT or Multiprocessing
355
356 Functions for calculating:
357 \int_a^b f(x) dx
358
359 Often not possible analytically, thus numerical aproach. Most differential
360 operations are possible analytically.
361
362 Riemann integral: Numerical approximation.
363 You divide the interval into subintervals (t_1, ..., t_n).
364 Riemann sum:
365 S_n = \sum_{i=1}^n f(x_i) \Delta x_i
366 Approaching the integral as the interval approaches 0.0. The geometric
367 interpretation corresponds to rectangles with height f(x_i) and width \Delta
368 x_i. The function passes trough the middle of the top side of the rectangles.
369 The integral exists also for discontinuous functions.
370
371 Rectangle method: Stable, does no assumptions
372 \int_a^b f(x) dx \approx \sum_{i=1}^n f(x_i) \Delta x_i
373 Trapezum method:
374 Instead of one point inside the interval, two points are calculated, and the
375 average is used.
376 \int_a^b f(x) dx \approx \sum_{i=1}^n (f(x_{i+1} )+ f(x_i))/2 * \Delta x_i
377 Adaptive methods:
378 The \Delta x_i is adjusted to how flat the function is in an interval.
379 As few points as possible are used.
380 Polynomial functions:
381 Analytical solutions exist. Taylor series expansion. w_i: Weight.
382 \sum_{i=1}^n w_i f(x_i) => \int_a^b f(x) dx
383 `---- quadrature -----`
384 n functions: {\phi_1(x),...,\phi_n(x)}
385 Often the polynomials are chosen: {1,x,x^2,...,x^{n-1}}
386 \sum_{i=1}^n w_i \phi_k(x_i) = I_k
387 I_k = \int_a^b \phi_k (x) dx
388 Gauss methods:
389 Also use the quadratures as tuning parameters, as the polynomial method.
390 w_i,x_i as tuning parameters: 2*n tuning parameters.
391 \int_a^b f(x) w(x) dx
392 Functions like 1/sqrt(x) or sqrt(x) are handled through the weight function.
393 See the lecture notes for details. Both nodes and weights are adjusted.
394
395 The method to use is dependent on the problem at hand.
396
397
398