Numerical Integration and Optimisation Theory

You can either open the above folders and read the descriptions inside or consult the 'Report.pdf' (it's very detailed, but you should somehow manage to read Russian) or just read the text below :)

//-------------------------------------------------------------------------------------------------------------------------

The first cluster of programs represents basics of numerical analysis:

1. The first piece of code, 'RecurrIntegr-Darboux.cpp', demonstrates instability occuring whilst calculating the integral $I_n = \int\limits_{0}^{1} \frac{x^n}{x + 6} dx = \frac{1}{n} - 6 I_{n-1}$ by a direct recurrent method, given that $I_0 = \log\left(\frac{7}{6}\right)$. Namely, calculating $I_{31}$ in such way gives us a cumulative error of order $\approx (-6)\cdot 10^{7}$, which in comparison to the real exact value of $I_{31} \approx 0+$ clearly persuades a courteous reader that this method makes no sense in our case. However, reverse method of calculation $\left( I_{n-1} = \frac{1}{6} \left( \frac{1}{n} - I_n \right) \right)$ quite expectedly turns out to be stable, and, given a value of $I_{62} \approx 0$, chosen for the sake of the number of iterations being the same as in the previous method, one can easily see that it works pretty much precise, yielding us $I_{31} \approx 0.00463$. The last part of this task was to compute the notorious $I_{31}$ by the definition of integral via lower and upper Darboux sums (https://en.wikipedia.org/wiki/Darboux_integral) $s_{31}$ and $S_{31}$, respectively. Of course, being unable to take the actual infinite limit, we will be dealing with the sequence elements close enough to be practically indistinguishable from the actual integral value. Uniform partition of $\left[ 0, 1\right]$ segment by $1000$ parts gives us $s_{31} \approx 0.0044$ and $S_{31} \approx 0.00456$, which is again pretty much what we anticipated it to be.

2. 'Analytic-Numeric-Diff.cpp' exhibits the non-trivial phenomenon of choosing $h$ which stands for the step of numerical approximation for differentiation, i.e. $\frac{f(x_0 + h) - f(x_0)}{h}$ instead of the actual analytical value of $f'(x_0)$ . For this purpose an elementary function $f(x) = 7 \sqrt{x^2 + x + 1} + 9 e^{-x} \sin x - \frac{1}{\log^2 x}$ has been chosen, with $x_0 = \pi$ and $f'(x_0) \approx 6.845544$ . The main part of the program constitutes an iterative loop with respect to $h$, on each step making $h$ lower by dividing it by $10$, starting from $h = 1$ and ending on $h = 10^{-20}$ . On each step the error $R_1 = \left| f'(x_0) - \frac{f(x_0 + h) - f(x_0)}{h} \right|$ had been calculated, put into a file and then plotted on the plane $h - R_1$ via gnuplot. Logarithmic scales were chosen for convenience, as the values vary much in magnitude.

3. 'Euler-method-1.cpp' is the first look at the brilliant numerical integration method known as Euler method (or rather first-order Runge-Kutta method), which, however, is ought to be treated accordingly to its peculiarities.

   a) The task was to solve numerically the following Cauchy problem: $\dot{y} = 2t, \quad y_0 = y(t_0) = 0, \quad t \in \left[ 0 , 10 \right].$ Three even partitions were used, with steps $h = 0.1, 0.01, 0.001$. Given the analytical solution $y(t) = t^2$, one can easily calculate the error for each step and see if the method works at all and where the error does maximise. Another gnuplot graph shows it clearly.

   b) The second part of the task was to solve another Cauchy problem: $y' = f'(x),\quad y(x_0) = f(x_0),\quad x \in \left[ x_0 , X \right] $, where $f$ is taken from the previous program, and $\left[ x_0 , X \right]$ is chosen so that the function is differentiable (and therefore integrable) on the entire segment. I picked $x_0 = 0.5$ and $X = 0.9$ since the function is discontionuous at $x = 1$ . This fact certainly has its effect on the quality of the numerical approximation, so we can see the error $R_1$ with the naked eye.

4. 'Euler-method-failure.cpp' represents the aforementioned peculiarities that should warn one from using it even for simple problems like solving classical oscillator ODE system: $\dot{x} = z,\quad \dot{z} = -x,\quad x(0)=0 , \quad z(0)=1$ . For proving this claim I obtained a table of the divergency values for numerical and analytical solution for various periods of time $T = 2\pi, 10\pi, \dots , 10^5 \pi$ and steps $h = 0.1, 0.01, 0.001$ . Furthermore, these tables were used for plotting the graphs of the actual values to have an obvious example of the stiff divergency. It was quite an amusement, I would say. Have a look and enjoy the power of Maths! :)