Numerical methods in C++

Markus Buchholz
Jun 21 · 7 min read

Numerical methods are typically used to solve mathematical models of nature and physical phenomenas. Each problem can be solved precisely. In this case we can say we solve the problem analytically or our solution are approximate. Previously the numerical method were used to solve very simple mathematical models, however increase of computational power (like FPGA, CPU, GPU or TPU) influenced increase of utilization of numerical computation to solve sophisticated models, described sometimes by high dimensional nonlinear equations.
When the computer(s) is (are) used to find a solution for a specific mathematical model (problem which can be modeled using mathematical equations) we can say that numerical method are used.
Numerical method in C++ will be discussed in two separate articles. Following part will cover the basics. The other part will cover more advanced aspects of numerical computation domain.

In following article I use the matplotlib library for c++ (clone repo). It is very simple library which allows you to plot data. The library is an API to very known Python library. The installation is not required, just only one header file.
In order to use the library, you need to first check your version of Python and path. You need to install also for the version of Python you have (for example using pip) the matplotlib.

#on terminal type:~$ whereis python   //you receive version of Python and paths:~$ pip3 install matplotlib //Python 3:~$ pip install matplotlib //Python 2

Each program you compile easily, but remember (in presented case) the header file has to be in the same folder as you cpp. You can modified the header path also.

//compile
g++ my_prog.cpp -o my_prog -I/usr/include/python3.8 -lpython3.8
// run
./my_prog
//folder tree.
├── my_prog
├── my_prog.cpp
├── matplotlibcpp.h

ITERATION

One of the most used operation used in numerical computation is iteration. The iteration, which you know exactly, can be presented as a set of the same operations repeated up the terminating conditions are satisfied. Frequently the conditions which we set are given by defining the precision (error) of the expected solution.

Iteration methodology in commonly used for finding the root of equation. Here I will present commonly use Newton’s method. The basic version of presented method can be defined as follows (the Latex version of all equation presented in following article you will find on my gist).

As you can see we find the next value of the function’s derivative f ′, and a previous root of f. We run the iteration process while the termination condition (precision) is satisfied.

Following example gives you opportunity to be familiar with discuss method. The example solves the simple quadratic equation.

Please refer to the example and adapt to you conditions.

We can also the plot how the program approaches the solution.

We can also apply to more complicated equation like high order polynomial (here third order):

We can also the plot how the program approaches the solution.

INTERPOLATION

Most function which are stored in numerical form (for example result of the computation stored on your PC) has discrete values for the specific domain of regarding function. It means the we have values of the function only in specific x(n) and x(n+1) points but no values between these two (for example x(0.5n).
Mathematical method to compute this values is called the interpolation. There is also extrapolation, which can be explained as an estimation— finding new values based on existing ones.

On below figure you can see discretize sin(x) function. The function is determined (has value) only in certain point. We will interpolate the other function and find the values between.

For interpolation we can apply mainly four different approaches: Newton, Lagrange’a, Chebyshev and spline interpolation.
Following example describes Lagrange’a method.

Following wikipedia we can find polynomial based on know values in nodes. We use following formula and compute n — order polynomial.

|i | x_i |f(x_i)|
|---------------|
|0 | 1 | 1 |
|1 | 2 | 3 |
|2 | 4 | 7 |
|3 | 8 | 11 |

For these certain (know) node values f(x_i) we can compute (using above Lagrange’a method) interpolated values. Let us assume, we want to know value for x = 6. We use the algorithm as follows:

f(6) = f(1)*[(6-2)(6-4)(6-8)]/[(1-2)(1-4)(1-8)]+
f(2)*[(6-1)(6-4)(6-8)]/[(2-1)(2-4)(2-8)]+
f(4)*[(6-1)(6-2)(6-8)]/[(4-1)(4-2)(4-8)]+
f(8)*[(6-1)(6-2)(6-4)]/[(6-1)(6-2)(6-4)] = 10.0476
// you can check computed value by running below code (values for different node will be computed and displayed)

The source code for discuss method is depicted below. Please, comment and uncomment the section for specific function (one presented in above example, sinus function and quadratic).

Below curves have been created by running last source cpp code. Please, note that blue curve represents know nodes, however the red one are computed — interpolated using Lagrange’a method). Presented results for some curves can not be very satisfactory, especially at the end of interpolation. The problem is related with applying interpolation high order of interpolating polynomials. More about this phenomenon is discuss here Runge’s phenomenon.

NUMERICAL INTEGRATION

Numerical integration is a set of method allowing for calculating the numerical value of define integral. Having a function f(x) we are able to estimate the value of this function in defined interval.
Most algorithms, which are used to compute numerically the value of integral use additives feature of the function, which can be written as follows:

As you can see we can split the defined function interval for small interval, compute integral (for shorter interval) and at the end apply a summation.
If following article we are going to use to different algorithms: rectangular method and trapezoidal method.

In rectangular method the function’s interval is divided by small rectangles. Our job is to for which we calculate just the area and at the end summ all computed rectangle areas. The discussed method uses following formula to compute the integral for certain function and interval.

More precise method applies not rectangles but trapezoids. In this can following numerical formula has to be applied.

Above equation has been applied in following code. The program computes the value of pi using the equation for area of circle with radius equals 1. I depicted two discussed methods, which can be plotted one by one. If so please comment and uncomment certain part of code in following file.

Below plots depict how the algorithms approached the solution with 9 numbers precision after comma.

NUMERICAL DIFFERENTIATION

Numerical differentiation is a set of methods to compute the value of derivative in defined point. Numerical computation in this case is willingly used, where there exist sort of difficulty to compute derivative for certain function. We deploy also numerical methods, where we know the values of the function only for discrete values of the function domain.
The easiest method to compute derivative is know as a Newton method which can be expressed as follows (in our case we can approximated derivative by computing the central difference quotient, h — is a small difference quotient):

Below code implements discussed Newton method and computes derivative for the polynomial third order. Please feel free to modify the function of your choice.

Thank you for reading.

Geek Culture

Proud to geek out. Follow to join our 1M monthly readers.