GopherProxy

	tnotes.rst - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	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 …
	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;…
	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…
	53 command, compiling the source code file with the same basename. `make -p…
	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 tabl…
	61 datapoints in an interval, the complete analytical function can be found…
	62 taylor-series of _all_ derivatives in the interval.
	63
	64 k-Polynomial: k+1 unknowns (variables). High polynomials tend to oscilla…
	65
	66 Cubic interpolation: The interval between points can be described by a f…
	67 of three polynomials.
	68
	69 Spline interpolation: Also made of polynomials. First order spline (a0+x…
	70 simply connects the data points linearly. The splines for each inter-dat…
	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 po…
	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, rememb…
	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 t…
	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}, …
	105 turn, c_{i+1}. Through recursion you find all unknowns by stepping forwa…
	106 The backwards solution can be found from the last data-point (y_n) by so…
	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…
	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 value…
	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 sim…
	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 …
	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 a…
	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 siz…
	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…
	245 sweep, which is faster than the classic rotation.
	246 Cyclic: Zero all elements above the diagonal, and do your rotations unti…
	247 matrix is diagonal. The matrix is converged if none of the eigenvalues h…
	248 changed more than machine precision. You will destroy the upper half plu…
	249 diagonal. If you store the diagonal in another vector, you can preserve …
	250 matrix values.
	251 In the beginning, you look for the largest element in each row, and crea…
	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 th…
	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 towa…
	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 …
	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…
	291 point. You sample your function, and fit a polynomial. Higher order poly…
	292 tend to oscillate.
	293 When solving ODE's the many sampled points create a tabulated function. …
	294 want to do something further with the function, you interpolate the poin…
	295 e.g. splines.
	296 Only three steps are needed if the equation is a parabola. Other functio…
	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 nee…
	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…
	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 …
	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 2c(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…
	324 e = \|\|y_1(full step) - y_1(2 half steps)\|\|
	325 You can also instead use RK1 and RK2, and evaluate the difference betw…
	326 two for the same step.
	327
	328 The effectiveness of the driver and stepper is determined by how many ti…
	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 t…
	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…
	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 differenti…
	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 \D…
	368 x_i. The function passes trough the middle of the top side of the rectan…
	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, a…
	375 average is used.
	376 \int_a^b f(x) dx \approx \sum_{i=1}^n (f(x_{i+1} )+ f(x_i))/2 * \De…
	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 metho…
	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 fun…
	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