C++ is a great programming language to calculate engineering problems. It is one of the best for those operations. Delphi is also a faster-compiled programming language that you can use in engineering problems. C++ compiler includes libraries designed to increase your productivity. We provide powerful and useful classes for strings, conversions, file IO, streaming, JSON, REST, and anything else you may require. In this post, we will explain how to solve millions of unknowns in Equations using the “Successive Over Relaxation Iteration” Method. In this post we will explain how to solve millions of unknowns in Equations by using the “Successive Over Relaxation Iteration” Method.

Table of Contents

## First, the science of unknowns in equations

Please take up a comfortable position and wear the appropriate safety equipment; we’re about to get all mathematical!

In mathematics there are Implicit Equations and we engineers mostly face solving numerical equations about the real life problems. Thermodynamics, Fluid Mechanics, Mechanism Theories, Structure Analysis, and many other engineering areas requires us to solve many *unknown* equations. Sometimes calculations these should be in real time or it should be simulating time fractions very precisely. Generally these implicit equations are exposed on a 1D, 2D or 3D grid examples (i.e compressed on uncompressed fluid flow, thermal distributions ) and these grids has nodes with parameters (i.e. Velocity, Pressure, Temperature etc.)

In mathematics, an **implicit equation** is a relation of the form *R*(*x*_{1},…, *x _{n}*) = 0, where

*R*is a function of several variables. In example, the implicit equation of the unit circle is is

*x*

^{2}+

*y*

^{2}− 1 = 0. More complex equations can be given as equations about to solve first and second law of thermodynamics, Momentum Equations, Naiver-Stokes Equations etc.

## Implicit functions and the Gauss-Seidel method

An **Implicit Function** is a function that is defined by an implicit equation, that relates one of the variables, considered as the value of the function, with the others considered as the arguments. These Implicit Equations can be solved by using Iteration Methods. One of the most popular equation is SOR Successive Over Relaxation Iteration Method (SOR). This method can be used to solve problems on 1D, 2D and 3D problems. In Numerical Linear Algebra, the **Successive Over-Relaxation (SOR) Method** is a variant of the Gauss–Seidel method for solving a linear system equations, resulting in faster convergence. A similar method can be used for any slowly converging iterative processes.

## One Dimensional Grid Problems: How do I solve a 3 Diagonal Matrix

Implicit Equations on 1D Grids produces Triangular Matrix forms. If the problem is on a 1D grid (for example heat transfer on a wire is 1D problem) this will produce 3 diagonal matrix forms. Because each 1D grid nodes on this wire connected to with left node and right except the first and the last one. That means there will be 3 variables (Left, Mid, Right) in each lines of the 2D matrix affecting to 3 Unknowns that we want to solve.

Simply equations in matrix form can be written as below.

In this equation;

A is 2D matrix (2 bars shows that it is 2D matrix),

U is 1D matrix form of unknown parameters from the each node,

q is 1D matrix form of right side of the equation.

Tridiagonal A matrix, U and q matrixes can be shown as below here,

**Two Dimensional Grid Problems**: How do I solve a 5 Diagonal Matrix?

If the problem is on a 2D grid (for example heat transfer on a plate is 2D problem) these equations will produce 5 diagonal matrix forms. Because each 2D grid nodes on this plane connected to with left, right, up and down nodes (or they are called East, West, North and South) except the corners and the edges. That means there will be maximum 5 variables in each lines of the 2D matrix.

In this equation;

A is 2D matrix (2 bars shows that it is 2D matrix),

U is 1D matrix form of unknown parameters from the each node,

q is 1D matrix form of right side of the equation.

5 diagonal A matrix, U and q matrixes can be shown as below here,

**Three Dimensional Grid Example**: How to solve a 7 Diagonal Matrix

As same examples above, if the problem is on a 3D grid (for example heat transfer on in a cube is 3D problem) this will produce 7 diagonal matrix forms. Because each 3D grid nodes on this cube connected to with left, right, front, back, up and down nodes (or they are called East, West, North, South, Up and Down) except the faces, corners and the edges. That means there will be maximum 7 variables in each lines of the 2D matrix.

**SOR Iteration Method **example of solving a 3 Diagonal Matrix

To understand this method, let’s go simple and we will use 3 diagonal matrix arrays. Create a new C++ Builder VCL Project and add two buttons (TButton), captioned ‘Create’ and ‘Solve’, and add a Memo (TMemo) to see outputs in a simple way.

Let’s define Number of Nodes, NN is 10 and let’s create our matrix arrays, m[][] is our 2D matrix (we said A before), u[] is for unknown values, x[] is the right side matrix (we said q before), and finally we define sol[] to hold solutions in iterations, so let’s define them in our VCL project,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
//--------------------------------------------------------------------------- #include <vcl.h> #pragma hdrstop #include "SOR_Iteration_Unit1b.h" //--------------------------------------------------------------------------- #pragma package(smart_init) #pragma resource "*.dfm" #define max(__a,__b) (((__a) > (__b)) ? (__a) : (__b)) #define min(__a,__b) (((__a) < (__b)) ? (__a) : (__b)) #define NN 10 // Number of Nodes TForm1 *Form1; double m[NN+1][NN+1], u[NN+1], x[NN+1], sol[NN+1]; //--------------------------------------------------------------------------- __fastcall TForm1::TForm1(TComponent* Owner) : TForm(Owner) { } |

Let’s assume that we know unknowns, solution is 1,2,3,4, and so on…, we can create a random equations with these solutions in a 3 diagonal forms as below,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
void generate_random_3diag_equations() { int i,j; srand( time(0)); for(i=0; i<NN; i++) u[i]=i; // First Line, Top Left Corner m[0][0] = (10+rand()%20)/100.0; // central node m[1][0] = (10+rand()%20)/100.0; // right node x[0] = u[0]*m[0][0] + u[1]*m[1][0]; // u0*E + u1*F for(j=1; j<NN-1; j++) { m[j-1][j] = (10.0+rand()%20)/100.0; // left node m[ j ][j] = (50.0+ rand()%10)/100.0; // central node m[j+1][j] = (10.0+rand()%20)/100.0; // right node x[j] = u[j-1]*m[j-1][j]+ u[j]*m[j][j]+u[j+1]*m[j+1][j]; } // Last Line, Bottom Right Corner //j = NN-1; m[j-1][j] = (10+rand()%50)/100.0; // left node m[j][j] = (10+rand()%50)/100.0; // central node x[j] = u[j-1]*m[j-1][j]+ u[j]*m[j][j]; } |

We can display all matrixes in this procedure by using `FloatToStrF() `

formating,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
void display_matrixes() { String str; for(int j=0; j<NN; j++) { str=""; for(int i=0; i<NN; i++) { str = str+ FloatToStrF( m[i][j], TFloatFormat::ffFixed, 9, 3 ) + " "; } Form1->Memo1->Lines->Add("|"+str+ "| |" + FloatToStrF( u[j], TFloatFormat::ffFixed, 9, 3 )+ "| |" + FloatToStrF( x[j], TFloatFormat::ffFixed, 9, 3 )+ "| |" + FloatToStrF( sol[j], TFloatFormat::ffFixed, 9, 3 )+ "|" ); } } |

## What is the Successive Over-relaxation Iteration method?

Successive over-relaxation iteration method is explained well here on Successive Over Relaxation Iteration Method. In this method there are some constants;**max iteration**: maximum number of iterations to exit from calculation, we took 1000**omega** : relaxation factor, we assumed 0.5 here**min_error **: minimum error allowed, we take 0.0001

In C++ Builder we can use this procedure below for the matrix forms,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 |
unsigned int calc_SOR( const unsigned int max_iteration ) { const double omega = 0.5; const double min_error = 0.0001; double err; double prev, next, psol; unsigned int iteration =0; do { err=0; for(int j=0; j<NN; j++) { prev=0; next=0; for(int i=0; i<j; i++) prev += m[i][j] * sol[i]; // these lines costly for(int i=j+1; i<NN; i++) next += m[i][j] * sol[i]; // these lines costly psol = sol[j]; sol[j]=(1.0-omega)*u[j] + omega*(x[j]-prev-next)/m[j][j]; err += pow( sol[j]-psol, 2); } err = err/NN; for(int i=0; i<NN; i++) u[i]=sol[i]; iteration++; Form1->Memo1->Lines->Add( UIntToStr(iteration) + ":" +FloatToStrF( err, TFloatFormat::ffFixed, 12, 6 ) ); }while(min_error<err && iteration<max_iteration); return iteration; } |

The example above uses full matrix and this is not wise for bigger memory allocations. We can use this matrix form below to calculate all diagonals for 1D, 2D and 3D problems.

1 2 3 4 5 6 7 8 |
struct matrix_line { float A,B,D,E,F,H,I; // values in a line of 7 diagonal matrix form to solve 1D, 2D and 3D problems short int a,b,d,f,h,i; // X positions (index) of A,B,D,E,F,H,I avlues }; struct matrix_line mx[NN+1]; |

But to understand this method first lets use full matrix forms, we may explain that struct based matrix form in the future posts. At the end, we can use these methods above in the `OnClick()`

event of our ‘Create’ and ‘Solve’ buttons. Let’s create our 3 diagonal matrix and display it in our Memo as below

## An example of a 3 diagonal matrix using C++

1 2 3 4 5 6 7 8 |
void __fastcall TForm1::Button1Click(TObject *Sender) { generate_random_3diag_equations(); Memo1->Lines->Add("Generated Equations in Matrix Form:"); display_matrixes(); } |

Let’s reset unknowns (we know the solution is 1,2,3,4 … so on) and let’s solve these equations by using SOR iteration methods.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
void __fastcall TForm1::Button2Click(TObject *Sender) { Memo1->Lines->Add("Reset Unknowns:"); for(int i=0; i<NN; i++) u[i]=0; display_matrixes(); Memo1->Lines->Add("Calculating..."); try { calc_SOR(1000); } catch(...) { }; Memo1->Lines->Add("Solution:"); display_matrixes(); } |

As a result, output will be as follows,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 |
Generated Equations in Matrix Form: |0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |0.000| |0.250| |0.000| |0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |1.000| |0.700| |0.000| |0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | |2.000| |2.200| |0.000| |0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | |3.000| |2.640| |0.000| |0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | |4.000| |3.810| |0.000| |0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | |5.000| |4.220| |0.000| |0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | |6.000| |6.020| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | |7.000| |7.180| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | |8.000| |6.370| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | |9.000| |3.740| |0.000| Reset Unknowns: |0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |0.000| |0.250| |0.000| |0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |0.000| |0.700| |0.000| |0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | |0.000| |2.200| |0.000| |0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | |0.000| |2.640| |0.000| |0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | |0.000| |3.810| |0.000| |0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | |0.000| |4.220| |0.000| |0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | |0.000| |6.020| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | |0.000| |7.180| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | |0.000| |6.370| |0.000| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | |0.000| |3.740| |0.000| Calculating... 1:14.927737 2:1.204483 3:0.110283 4:0.014905 5:0.003814 6:0.001554 7:0.000785 8:0.000444 9:0.000270 10:0.000172 11:0.000114 12:0.000077 Solution: |0.210 0.250 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |0.094| |0.250| |0.094| |0.160 0.500 0.100 0.000 0.000 0.000 0.000 0.000 0.000 0.000 | |0.944| |0.700| |0.944| |0.000 0.280 0.570 0.260 0.000 0.000 0.000 0.000 0.000 0.000 | |2.067| |2.200| |2.067| |0.000 0.000 0.240 0.560 0.120 0.000 0.000 0.000 0.000 0.000 | |2.953| |2.640| |2.953| |0.000 0.000 0.000 0.220 0.500 0.230 0.000 0.000 0.000 0.000 | |4.031| |3.810| |4.031| |0.000 0.000 0.000 0.000 0.140 0.540 0.160 0.000 0.000 0.000 | |4.996| |4.220| |4.996| |0.000 0.000 0.000 0.000 0.000 0.240 0.570 0.200 0.000 0.000 | |5.988| |6.020| |5.988| |0.000 0.000 0.000 0.000 0.000 0.000 0.230 0.520 0.270 0.000 | |7.018| |7.180| |7.018| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.100 0.540 0.150 | |7.990| |6.370| |7.990| |0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.220 0.220 | |9.013| |3.740| |9.013| |

As you see in this example we achieved to results about 0, 1, 2,3, 4.. in 12 iterations with the error 0.000077. You can use this method to solve millions of unknowns on your grid systems, i.e on fluid mechanics, thermodynamics, stress analysis, and to solve many other linear equation problems.

## More about iteration methods

Iteration methods are good to solve some problems, also sometimes hard to understand and hard solve some problems because of their behavior moves from convergence iterations to divergence iterations. That causes wrong solutions. So knowing these methods is not enough, also you should know well about it’s behavior of relaxation factor, error factors, assumptions and values in matrix also important. Sometimes you must use dimensionless parameters, or reformed matrixes to achieve results or to have lower iterations.

In some random generation tests, you will see that errors in iterations goes divergence instead of being convergence. There are SOR Revised ad some other new methods to solve these problems. Please check academic papers and latest releases for better solutions. In that point you will see that knowing mathematics well and knowing programming well is not enough to solve these kind of problems. You must also know the behavior of these kinds of methods.

**RAD Studio C++ Builder is a great environment for learning to use C++ and is also powerful enough for all your development needs. Why not download and try C++ builder today?**

Design. Code. Compile. Deploy.

Start Free Trial Upgrade Today

Free Delphi Community Edition Free C++Builder Community Edition