Complex numbers and equation optimization using GA

[Nighteyes]

Hi
I am currently trying to write a genetic algorithm to find the zero off a very big equation.

The problem is, is that it involves numerus complex numbers. Only at the very last stage in the program do I want to convert these complex numbers into real numbers.

This is done so I can plot the real and imaginary parts of the eigenvalue pairs so I can find the point were they cross the axis.

BUT I can not find away of of using complex numbers in c++. and to convert these equations into two X = Real(X) + i * imaginary(X). would literally take days of rearranging equations !!.

Is there any way I can get around this problem?

Also whilst im at it my idea was to then import the complex numbers Eign_1_1 -> Eign_2_2 into another program were I would then split them into Real(Eign_1_1) and imag(Eign_1_1) which the zero would be then found using an evolutionary algorithm.

Any ideas on how to do this part?

Thanks

Craig

[Philip Nicoletti]

Simple example:

#include <complex>   // not complex.h
#include <iostream>

using namespace std;



int main()
{
    complex<double> a(1,1);
    complex<double> b(3.2,-3.4);

    complex<double> c = a + b;
    complex<double> d = sqrt(a);

    cout << a << "\n";
    cout << b << "\n";
    cout << c << "\n";
    cout << d << "\n";

    return 0;
}

[dude_1967]

Yes, like Philip said, you must become adept with the use of the STL template class called complex, located in <complex>. You will note that the implementation of complex is well thought-through. It offers simple member functions real(...) and imag(...) for extracting the real and imaginary parts.

In addition, if I may ask, what is the implementation of the Bessel functions? Are you using GCC or the implementation in cephes or something like that?

I feel as though you are getting off to a less-than-optimal start with your C++ design. In your sample code there are several data tables and constants which could be wrapped up nicely in STL containers.

Furthermore, there are many repeated code sequences and repeated calculations---probably very poor in run-time as well as challenging for an optimizing compiler, in general, though, strongly ill-advised.

I could offer clear suggestions for a better design of your numerical architectures if you would like to take a bit more time and learn any of these. Unfortunately my time is rather tight, but it would be a shame to let another numerical algorithm developer get away with writing out-dated FORTRAN77 styles within a weak C++ design concept. I do not wish to be mean or cruel with these judgements, but rather alert you to the fact that things can be done a lot better, more stable and elegant, run-time optimized, etc with the right programming idioms.

Sincerely, Chris.

[Nighteyes]

DUDE!! ...1967 Yes please give me all the critisim and advice possible please.

I am completely new to C++, got "object orientated programming in C++" last firday.

And yeah that was origonally writen in matlab. I treid doing it properly in c++ the first time I tried to programm it (yesterday) but could not work out how to do it better. (I tried to do it all in a class..... needless to say it didnt work)


// useless background below
At the moment my undergraduate project is on this code which is based on a paper by Holt Ashley called "Role of Shocks in the "Sub-Transonic" Flutter Phenomenon" but there are alot of errors in the paper so i cant get it working.
So instead of spending another year trying to get it working I though I might be able to get my prof. to allow me to change the project to do with something about optimizing using evolutionary programing.
As the method in the above program is the same as that used by uni called the BIFOR code I figured that if i can get the genetic algrithim to converge Eign_X_X to zero faster than the current mehtod he would allow me to change my project.
Also I have a parametic C++ code for an aerofoil which gain I think could be optimized using evolutionary programming.

Hence im learning C++.

So any and all help anyone can give me to make my programs run fast/ be more efficent will be gratefully resived.

OH and the bessel functions, yeah... Well as I mentioned in the background section the equations are all form the Holt Ashley paper. The F1, F2 and F3 terms are part of his linear aerodynamic solition. So I am unable to get rid of them. And the hn function is for theodorsans function (equation called theo) can i cant change that either unless I rewrite the rule book for aeroelastic problems

So in short there only in this code cause some proffesor decided they would be good :S

Thanks

Craig

[dude_1967]

OK Craig,

It is late in Germany now. So I'll only be able to dive into this tomorrow morning. I will look at your sample file once again and try to work out some sensible suggestions.

Could you in the meantime try to answer these questions?

- What compiler system are you planing to use?

- Where are the Bessel functions defined? For example where is the subroutine j0(...) actually implemented?

- Are you planing to use spherical Bessel functions of half-order (sometimes known as Ricatti Bessel functions) or normal Bessel functions of integer order?

- Please try to describe the planned use of complex numbers. I have looked at your code and seen nothing but reals.

- I am not familiar with your Eigenvalue method. Do you feel that you will need any kind of root-finding method such as the secant method, or is the root-finding the point of the Eigenvalue algorithm.

- Are you familiar with this reference: http://www.amazon.com/Object-oriente...6327189&sr=1-1
It should be available at one of your libraries.

Sincerely, Chris.

[Nighteyes]

ok i'm at uni at the moment so cant give you good answers till tonight.

1) using G++
2) the bessel functions coming straight from cmath.h.
3) normal Bessel functions of integer order

the rest will have to wait till i can give it a good look at an try and work out the answers.

also I will try and get hold of that book.

Thanks

Craig

[dude_1967]

OK Craig,

Based on your answers I took another look at the sample code.

We need to go through a few points.

Point 1:

First of all you need to gain an understanding of <complex< numbers in C++. There is a template class called std::complex which implements the mathematics of complex numbers. The use is quite simple:

#include <iostream>
#include <complex>

int main(void)
{
  std::complex<double> a(1.1, 2.2);
  std::complex<double> b(0.5, 0.6);

  std::complex<double> c = a + b;

  std::cout << c.real() << std::endl;
  std::cout << c.imag() << std::endl;
}

So when you take the square root of the double -1 in your program this will not work at all. You seem to be writing things in code like x + i*y. This is not at all correct. Rather, you define a complex number and use normal mathematical expressions on it. You need to familiarize yourself with the complex number template in the C++ STL.

Point 2:

I always establish a sensible naming convention for the special functions of pure and applied mathematics such as the Bessel and Hankel Functions that you need to use. I always use the naming convention and parameter organization of Mathematica. I find the 'usual' symbols j0 or y0 which are commonly used for Bessel functions in Posix to be non-intuitive and confusing. So I always wrap these in a function like the ones shown here:

#include <cmath>
#include <complex>

double BesselJ(const int n, const double x)
{
  if(n == 0)      { return ::j0(x); }
  else if(n == 1) { return ::j1(x); }
  else            { return ::jn(n, x); }
}

double BesselY(const int n, const double x)
{
  if(n == 0)      { return ::y0(x); }
  else if(n == 1) { return ::y1(x); }
  else            { return ::yn(n, x); }
}

std::complex<double> HankelHl(const int n, const double x)
{
  return std::complex<double>(BesselJ(n, x), BesselY(n, x));
}

It is interesting to note that the Hankel function of the first kind actually returns a complex type.

Point 3:

The next thing that I noticed in your code is that you force the compiler to repeat many calculations. This is poor design and results in very poor run-time characteristics. For example you compute the factor 1-M*M and its square root very often. In addition you repeatedly call j0(lambda) when, in fact, one call would be sufficient.

Point 4:

Try to use sensible text allignment schemes to improve the readibility of your code.

Point 5:

We usually discourage the use of 'using std' because this can generally lead to confusion with naming. We suggest that you use std::complex, std::cout, etc.

Point 6:

We need to take a good look at the kind of complex mathematics that you really want to do. By looking at your calculation of the variable 'theo' it seems like you might actually need to compute Bessel functions of complex argument. This is not available from g++ directly. You would need to find an implementation for this. I could help you with this. But first we need to review the complex mathematics which you need to use. Try to describe this and, if this is not sufficient, then we actually need to look at the equations in the literature.

Point 7:

If you are going to work with me then you are going to use spaces and no tabs. I use two spaces for indentation.

----------------------------------------------------------------------------------

So I put all this together and wrote a suggestion on how you might want to start with this kind of thing. It is based on your original code but only covers about 1/4 of it.

You need to study these things.

Take a look at my suggestions. Try to answer my previous questions as well as the questions in this post. And then we can move on.

Below is a code sample based on your original work.

Sincerely, Chris.

#if defined(_MSC_VER)
  #pragma warning(disable:4996) // Disable Microsoft compiler warning
#endif

#include <iostream>
#include <cmath>
#include <complex>

double BesselJ(const int n, const double x)
{
  if(n == 0)      { return ::j0(x); }
  else if(n == 1) { return ::j1(x); }
  else            { return ::jn(n, x); }
}

double BesselY(const int n, const double x)
{
  if(n == 0)      { return ::y0(x); }
  else if(n == 1) { return ::y1(x); }
  else            { return ::yn(n, x); }
}

std::complex<double> HankelHl(const int n, const double x)
{
  return std::complex<double>(BesselJ(n, x), BesselY(n, x));
}

int main(void)
{
  static const double k = 0.5;

  ////////////////// initialization /////////////////
  const double M     =  0.90;  // Flow field mach number
  const double rho   =  0.50;  // ambient air density
  const double Xc    =  0.00;  // dimensionless location of cg
  const double Xs    =  0.00;  // dimen. location of shock foot	
  const double a     = -0.20;  // dimen. location of elastic axis
  const double b     =  0.75;  // half cord 
  const double omega =  2.00;  // omega w / omega alpha
  const double mu    = 50.00;  // mass ratio

  ////////////////////////////////////////////////////
  ////////////// variables of F functions ////////////
  const double m2                = M * M;
  const double one_minus_m2      = 1.0 - m2;
  const double sqrt_one_minus_m2 = ::sqrt(one_minus_m2);

  const double kdash  = k / one_minus_m2;
  const double lambda = kdash * (m2 - (sqrt_one_minus_m2 * ::log(1.0 + (sqrt_one_minus_m2 / M))) - log(M / 2));

  const std::complex<double> hankel_one_kdash = HankelHl(1, kdash);

  // TBD: Not yet finished.
}

[Nighteyes]

wow! Thanks mate, cant wait till I get home will have a good look though that!

I know the answers to the other two questions in the previous post, but it will require a proper explanation (more than I can give here with out my report at hand).

So will either re-compile my report to just include the relivent bits or pm you the pdf and refer you to the relivent page + an explanation put into context of c++

Probably wont be for another 4 hours though.

Thanks

Craig

[Nighteyes]

Right, the Eigenvalues.

This is a flutter problem solver... The equations eventually come down to the form

AA(omega_theta/omega)^4 + BB(omega_theta/omega)^2+CC = 0 (when flutter occurs)

we will refer to (omega_theta/omega)^2 as sigma form now on, so we have;

AA*sigma^2+BB*sigma+CC=0

were AA, BB and CC are some functions of k and are complex numbers

we then follow this procedure to solve the equations:

1) Specify a set of trial values for k values, usually from 0.001 to 1.0.
2) For each value of k calculate the aerodynamic coefficients.
3) Solve the flutter determinate, which is a quadratic equation with complex coefficients, for the values of sigma that correspond to each of the selected values of k. Note that these will be complex in general, the real part being an approximation of sigma, and the imaginary part being related to the damping of the mode.
4) Interpolate to find the value of k at which the imaginary part of one of the roots becomes zero. This can be done approximately by plotting the imaginary part s of both roots against k, so that the value of k at which one of the imaginary parts crosses the zero axis can be determined. This value of k then has a corresponding real value of sigma, which provides the value of omega.
5) Determine flutter speed U=(b*omega)/k
6) Repeat for different values of mass ratio, Mach number or frequency ratio to map out the flutter boundary with respect to one of the above.

See the attachment to for how the eigenvalue pairs look.

The use of the evolutionary algorithm will be to find the zero closest to k=0 on one of the eignvaule pairs.

So as you can see the use of these Eigenvalues is unavoidable.

Hope this helps.

Just finished going though your suggestions and have come across a few problems.

Firstly as I understand it to enter a complex number into c++ it requires separately entering both the real and the imaginary parts of a complex number.

If this is the case how do you multiple something by i in the middle of an equation?

I have also programed the rest of this up trying to prevent using divide and instead creating variables of say one_over_X = 1 / X; so it is possible to just multiple bu one_over_x in future limiting the uses of divide.

Is this a good thing to do?

Finally when i enter " std::complex<double> F1, F2, F3; " the compiler goes haywire??

well thats about it for now. Hope you think this is alot better than attempt 1... and no tabs used :P

Thanks

Craig

[dude_1967]

Craig,

I just finished going through your second round code.

We are coming along pretty well. But this will take a few more rounds.

We need to:
1) Fully understand the use of complex
2) Get the numerics correct
3) Get the right architecture

Let's concentrate on numbers one and two for now.

I have reworked your code and it does compile now using Visual C 9 as well as g++ 4.3.0.

To answer your question about the imaginary number i, you can simply define i to be (0.0, 1.0). I have done this in the reworked code and actually created a subroutine called I(...) which returns a constant reference to this imaginary number. Again, the name 'I()' is taken from the Mathematica naming convention. So this number can simply be used anywhere any complex can be used.

You had significant troubles because you were trying to set doubles equal to the values of complex'es. You must find the right place in the algorithms to switch over to real values by extracting the real and imag parts using the member functions real(...) and imag(...). In my re-worked code, the calculations are carried out to the end using complex'es.

Where do you want to switch to real values?

I am concerned with the complexity of the expressions for F1, F2 and F3. The lines which set F1, F2 and F3 are overly complicated. These lines should be simplified, you know broken down into multiple lines featuring their their constituent parts using additional helper variables, combining them later into F1, F2 and F3.

In general your code could be improved by making more use of parentheses and groupings. You make heavy reliance on operator precedence. This was especially evident for AA, BB, etc. I did a bit of grouping but, in doing so, probably actually broke the code.

You must go through all of the mathematics and we need to make sure that we are actually getting the right numbers.

Usually I start small, check each algorithm line, and progress to the final answer. We are a bit backwards here because a body of unchecked mathematical code exists and we are refactoring it. We still need to establish the numerical integrity top-down up to the end.

So I want you to do the following:

1) Take my re-worked code.
2) Simplify the things that I suggested and establish better use of groupings and parentheses.
3) Identify the right place to extract the real and imag parts.
4) Try to verify the numerical integrity of some intermediate results such as F1, F2, F3.

We will go a few rounds after that. I will be out and about a lot because it is the weekend. Perhaps I might not respond until Sunday. But I will respond.

Good luck with the next steps.

Sincerely, Chris.

[Nighteyes]

Thanks for the help again Chris.

I finally see how complex numbers works in C++, its so simple now

I had originally made a matlab program for this problem so it was easy enough to check all the numerics were correct. Though I have noticed that as the programe progresses rounding errors seem to creep in. I have attached a .csv file which shows the errors from Matlab to C++.

As can be seen most of the errors are less than e-07 with the largest being e-03.
This does not really matter to me as I going to be demonstrating the difference in convergence speed between to different methods. But how many decimal places does C++ carry forward?

Finally I will be converting the Eig_1_1 and Eig_2_1 terms into real and imaginary bits for use in the solver (going to be doing both the Newton-Raphson method and a Evolutionary algorithm). So were they are converted can either be at the end of this program or when they are imported into solver, which ever is easier/better.

Craig

p.s. had to rename Errors.csv to errors.txt. It is still comma separated though

[dude_1967]

I finally see how complex numbers works in C++, its so simple now

Yes, it is simple and powerful.

I had originally made a matlab program for this problem so it was easy enough to check all the numerics were correct. Though I have noticed that as the programe progresses rounding errors seem to creep in. I have attached a .csv file which shows the errors from Matlab to C++.

See my changes to the code using precision. Let's make sure that we are printing out all of the precision and then find out if there is a round-off error and, if so, then where does it come into the calculations. I compiled with two compilers and got very similar results to about 13 decimal digits. Please remember that it might be Matlab with the precision limitation. You can set the default precision of Matlab, but I forgot how to do this. You might also need to use quoted numbers in Matlab in order to maintain precision beyond the default printed precision which, I believe, might be only 6 digits.

But how many decimal places does C++ carry forward?

See my changes to the code using <limits>. The decimal precision of a C++ double is specific to the compiler and plattform. However, most common systems normally support 15 decimal digits of precision for a C++ double.

Craig, I see you are moving right along. So let's look at a few next steps.

I did the following (see my comments):
1) Introduced the optional use of C++ namespaces for architectural organization.
2) Set the default precision of std::cout to the decimal precision of double using the member function precision in combination with <limits>. This technique is valid ofr any kind of output stream including output files.
3) Set up a subroutine architecture and introduced the use of a C++ STL std::vector of complex'es. Take a look at the loop in main.
4) There were some unused variables. Compile with g++ -Wall -O3 test.cpp -o test.exe. The option Wall includes all warnings and showns, among other things, unused variables. The option O3 using the letter big-O stands for the level 3 optimization which will cut much faster code.
5) I used formated output and extracted the real and imag parts of the Eigenvalues. Take a look at the loop in main.

Take a look at these new things.
- Are the numbers are getting any better?
- Do you want to do even further investigations or use higher precision?

Sincerely, Chris.

[Nighteyes]

Thanks alot mate.

Yeah I looked at the numbers, they look fine. As the precision difference is much lower than the number of decimal places I will be using it does not matter. Also to find out exactly were the difference in precision appears will be to much of a detour to make it worth while at this time.

Though once I get my evolutionary algorithm working I definitely will be interested in methods to reduce this kind of error. Though my understanding of this is limited to something called truncation error. The best way I know of to reduce this is to multiple all numbers by say 1e12 to make numbers none floats then dividing by 1e12 at the end. (apparently this also speeds up cpu time (allegedly))

I am a little confused about "architectural organization" does this refer to computer architecture or is this about the architecture of the source code?

What is the effect/benefit of introducing this?

Also what is the benefit of introducing Eigen as a vector?

Reading though the code I think what it does is this
vector:
[Eigenvalue one
Eigenvalue two]

So the function can return more than one value. Is this correct?

I am now trying to export these Eigen values into a solver. So I wrote a Newton-Raphson method solver for this program (see attached). I have tried to incorporate the things you have taught me about programming into this program but there is little room for variance int this program due to its simplicity.

My only real concern is calculating the modulus efficiently, at the moment I do this by: sqrt((x*x)) which my not be the best method.

Could you please show me how to get the Eigenvalues into this program. I think I need to make the Ashley program into a h file, but my brother has my c++ book at the moment so Im stumped :S.

Thanks again

Craig

// Implimentation of Newton-Raphsons method
#include<iostream> // for cout
#include<cmath>  //for sqrt

double Ashley_Eig_1_imag (double k) // The imaginary part of Eignenvalue one from the Ashley program
{
double ans = k*k*k;  // this was just to test it as I do not know how to import Egin_1_imag
return ans;
}
int main()
{
  double fdash, xn, xn_pluse_1, mod_xn, mod_xn_plus_1_minus_xn; // variables
  int N=0;
  const double h = 0.000001;  // Not sure but I think the smaller h is the more accurate fdash is
  const double tol = 0.00001;  // tollerance for program ending
  std::cout << "Input Initial guess: ";
  std::cin >> xn;  // initial guess is important for N-R method

do
{
  fdash = ( Ashley_Eig_1_imag(xn+h) - Ashley_Eig_1_imag(xn) ) /  h; // differntation of Eigenvalue 1 imaginary part
 

  if(fdash == 0) // if this is true method will not work
  break;

  xn_pluse_1 = xn - (Ashley_Eig_1_imag(xn) / fdash);  // calculate true intersept point
  mod_xn = sqrt((xn*xn));  // calcuate modulus of xn (not sure if this is the best way to do this
  mod_xn_plus_1_minus_xn = sqrt(  ((xn_pluse_1 - xn)*(xn_pluse_1 - xn))  );


  if(mod_xn_plus_1_minus_xn <= (tol*mod_xn))  // end conition
  break;

  std::cout <<"("<< N <<")" << "Break func: " << mod_xn_plus_1_minus_xn << "<=" << (tol*mod_xn) << std::endl;
  xn = xn_pluse_1;  // set the new xn value as the old one reasy for next run though
  N++;
} while(N < 100000);  // max number of runs 100000
std::cout <<"("<< N <<")" << " K is: " <<  xn_pluse_1 << " at Zero" << std::endl;   // print out answer
return 0;  // cause all the examples of c++ code has this :D
}

[dude_1967]

Good. It seems like you are moving right along, Craig.

As the precision difference is much lower than the number of decimal places I will be using it does not matter.

How many decimal digits of precision do you need?

Though once I get my evolutionary algorithm working I definitely will be interested in methods to reduce this kind of error.

This can be a rich topic. We can get into it later, if need be.

I am a little confused about "architectural organization" does this refer to computer architecture or is this about the architecture of the source code?

Yes, I am talking about the organization of the source code.

What is the effect/benefit of introducing this?

Organizing program components into individual files, classes and namespaces can generally improve program readability and maintenance. In addition, it is the foundation for design. Your program architecture should represent, in some way, the structure of the problem at hand. Subroutines are created with sensible names and corresponding functionalities. These are organized in modules and namespaces, etc. Again, this is a rich topic about which many books can and have been written. You just have to work with these things, study available literature and modern programming techniques and thereby develop a good sense for program architecture and design.

Also what is the benefit of introducing Eigen as a vector?

There could have been many different solutions for this part. The point was to return the pair of complex values in some sort of single unit. I just happened to select std::vector<complex<double> > for this purpose. It could have just as well been C++ structure or a C++ class which embodies the two complex values. Vector is one of several container classes implemented in the C++ STL. There are also deque, list, stack, map, etc. Each container has its special advantages, i.e. speed, push-optimization, sort, etc. Vector and deque are often well suited for storing numeric results. There are good FAQ's, reports and a host of good literature describing STL, vector, etc.

So the function can return more than one value. Is this correct?

In this case, GA::Eigen returns a single vector containing two complex values. We also note here that returning a vector is a design choice which is poor for large vectors, but maybe OK for small vectors. For large data amounts, an input non-const reference to a vector would be better.

I am now trying to export these Eigen values into a solver. So I wrote a Newton-Raphson method solver for this program (see attached). I have tried to incorporate the things you have taught me about programming into this program but there is little room for variance int this program due to its simplicity.

Actually there is a lot of room for improvement in your Newton-Raphson concept. There is a very elegant implementation if Newton iteration in the celebrated Numerical Recipes in C++ book somewhere around Page 360. Newton iteration requires the evaluation of a function as well as its derivative at a point. I do not really like the organization of your Newton iteration. I would write it in a generic fashion using a generic function which calculates the value of a function and its derivative. Furthermore, I find that you should really use an improved method for calculating the derivative. You are using a first order rule with a very large step-size. This will not be very accurate.

My only real concern is calculating the modulus efficiently, at the moment I do this by: sqrt((x*x)) which my not be the best method.

Use ::fabs(x) for x a real double.

Could you please show me how to get the Eigenvalues into this program. I think I need to make the Ashley program into a h file, but my brother has my c++ book at the moment so Im stumped.

You need several files. Make files Ashley.cpp and Ashley.h. Put the function prototype for GA::Eigen in Ashley.h. Be sure to include the namespace and remove the static qualifier. Make another file called Newton.cpp and include Ashley.h in Newton.cpp.

1) You need to work on the Newton iteration. It is not very good. Go look at the implementation in "Numerical Recipes in C++" or a similar book in your library. Use a higher order method (maybe third or fourth order) for the derivatives combined with a step-size of, say 10^-3 or 10^-4.
2) Try to organize your program into several files and get back to us if you need additional help with this. You can build your program with a single g++ command line if you basically use the same command-line as we had and instead of using one cpp file, just list the cpp files.

Keep going. You are coming along well. These methods in C++ can be very effective for creating numeric software of high quality. And this will definitely aid your studies and career.

Sincerely, Chris.

[Nighteyes]

Thanks again Chris.

Guess it was just wishful thinking that the N-R program was good...

I am still waiting for the numerical methods c++ book, but I do have the C one which should do.

It will probably take me a while to have a really good shot at this (hopefully done by the end of the weekend) but will post back as soon as I am done.

oh and precision wise I think I need at most 1e-8 but can probably get away with 1e-6.

Craig