Problems 8

Theory.
1. How do you allocate GSL vectors and matrices?
2. Do you need to manually free the memory allocated for your vectors and matrices?
3. Which GSL functions can be used to add two gsl-vectors? Hint: there is more than one.
4. Which GSL functions can be used to multiply a matrix and a vector?
5. Which GSL functions can be used to add two matrices?
6. Which GSL functions can be used to multiply two matrices?
7. How does one solve a system of linear algebraic equations 'Ax=y' with GSL? Name as many options as possible.
8. What does the function gsl_blas_dgemm do?
9. After you call the GSL's LU decomposition routine on a matrix, is the matrix preserved or destroyed?
10. How can you find out the declaration of struct gsl_odeiv2_system ?
  Hints:
  - #include<gsl/gsl_odeiv2.h> and gcc -E;
  - gsl-config --cflags and view.
Practice.
1. 1. Solve, using GSL, the linear system
    
    [ 6 -5 -9 ] [x0] [-9]
    
    [ 2 -2 -4 ] [x1] = [0]
    
    [ 0 -5 4 ] [x2] [-8]
    
    Hint: you can solve linear systems using one of the following methods,
    - Householder solver for linear systems gsl_linalg_HH_solve
    - LU-decomposition using the GSL functions gsl_linalg_LU_decomp gsl_linalg_LU_solve
    - QR-decomposition using the GSL functions gsl_linalg_QR_decomp gsl_linalg_QR_solve
    See
    GSL manual → Linear Algebra → Householder solver for linear systems;
    GSL manual → Linear Algebra → LU Decomposition;
    GSL manual → Linear Algebra → QR Decomposition;
    GSL manual → Linear Algebra → Linear Algebra Examples.
  2. Check that you've got the correct solution.
2. Consider the equation of equatorial motion of a planet around a star in General Relativity,
```
u(φ)'' + u(φ) = 1 + εu(φ)² ,
```
  where u(φ) ≡ 1/r(φ), r is the (reduced-circumference) radial coordinate, φ is the azimuthal angle, ε is the relativistic correction (on the order of the star's Schwarzschild radius divided by the radius of the planet's orbit), and primes denote the derivative with respect to φ.
  1. Integrate this equation with ε=0 and initial conditions u(0)=1, u'(0)=0 : this should give a Newtonian circular motion.
  2. Integrate this equation with ε=0 and initial conditions u(0)=1, u'(0)≈-0.5 : this should give a Newtonian elliptical motion. Hint: u'(0) shouldn't bee too large or you will lose your planet.
  3. Integrate this equation with ε≈0.01 and initial conditions u(0)=1, u'(0)≈-0.5 : this should illustrate the relativistic precession of a planetary orbit.
  Hints:
  - You can rewrite the equation as a system of ordinary differential equations (ODE) by introducing the functions y₀=u and y₁=u', which gives the following ODE,
```
y₀' = y₁
y₁' = 1-y₀+ε*y₀*y₀
```
  - If you have calculated the "data"-file as two columns with φ and u, you can plot the orbit with the following gnuplot/pyxplot command,
```
 
plot "data" using (1/$2)*sin($1):(1/$2)*cos($1) with lines notitle 
```

[	6	-5	-9	]	[x0]		[-9]
[	2	-2	-4	]	[x1]	=	[0]
[	0	-5	4	]	[x2]		[-8]