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-oriented-Numeric-Computing-Scientists-Engineers/dp/0387989900/ref=sr_1_1?ie=UTF8&s=books&qid=1216327189&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:D

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 nowYes, 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. :thumb:

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

[dude_1967]

Guess it was just wishful thinking that the N-R program was good...Please don't be discouraged. In fact, I highly approve of your efforts and ability to learn.

When I was in grad school years ago, the general programming skills of university reserch professionals and professors were very poor. There has been unfortunately little improvement in this area due to the lack of communication between cutting-edge professional programmers and university researchers. We are trying to close this gap. And you should keep trying to improve as well.

I will be away from all offices on Friday and return on Sunday.

My last post was a bit terse on the description of the Newton iteration. If I get a chance, then I might be able to work out a sample program for computing the derivative of a generic function. This would provide you with a better impression of what I am talking about in terms of style and implementation.

Keep working on this stuff. You are going in the right direction. Post your progress and further questions when you are finished with the next steps in a few days.

Sincerely, Chris.

[dude_1967]

Hi Craig,

I was thinking about my recent response about Newton Iteration. At first you might want to stick with your first-order estimate of the derivative. This will keep the implementation simple.

That Numeric Recipes book that I had mentioned really does a good job explaining this kind of stuff.

If I get a chance, then I might be able to work out a sample program for computing the derivative of a generic function. This would provide you with a better impression of what I am talking about in terms of style and implementation.It was so much fun that I looked at some of my notes from 20 years ago. I wrote down the central difference rules for numerical derivatives and worked out a sample C++ file for testing these. They are really neat little formulae. See the attatchment.

Sincerely, Chris.

[Nighteyes]

wow thanks mate!

Got my C++ programming book back of my bro today so gonna read a few more chapters in that tonight, esp multi-file programs.

I am gonna stick to the first order calculation, just because this is what is used in the code I am trying to replicate to show genetic algorithms/ differential evolution algorithm can do it faster. (this may not be the case)

I am set on learning C++ to the best of my abilities and as fast as possible so I will not be discouraged, just encouraged to try an do even better with the next draft of the code

The "Guess it was just wishful thinking that the N-R program was good..." is my everlasting optimism that I can do something perfect first time (like that will ever happen :S)

I had a good look though your code for the derivatives and then looking at the book called "Numerical recipes in C" I have decided to change my fdash function from

fdash= F(x+h)-F(x) / h to fdash=F(x+h)-F(x-h) / 2h

So I can calculate the actual difference this makes on the calculation of the zero point although this will have to kept as fdash= F(x+h)-F(x) / h for calculating weather a Genetic algorithm or newton method is faster for getting convergence

also interestingly the book suggests by setting a new variable temp as

temp = x+h so temp is an exact representable number in binary

The error can be reduced significantly. Would setting temp, then for example 1, count as the smallest possible exactly represented binary number.

Also in you notes what is the symbol in 1st derivatives

fdash = m1 + unknown_symbol (dx^2)

Thanks again for all the help, I think im learning much faster than if I was doing this on my own.

Also how do you access g++ documentation/help files? cant find it anywhere on my computer, but then again I only started using linux and programming 2 weeks ago now

Craig

P.S do you think it is a good Idea to change the title of this thread to "Numerical programing" or something like that?

[dude_1967]

I am gonna stick to the first order calculation, just because this is what is used in the code I am trying to replicate to show genetic algorithms/ differential evolution algorithm can do it faster. (this may not be the case)That sounds fine.

I am set on learning C++ to the best of my abilities and as fast as possible so I will not be discouraged, just encouraged...Great! C++ is a great language for numerical calculations.

I had a good look though your code for the derivatives and then looking at the book called "Numerical recipes in C" I have decided to change my fdash function from

fdash= F(x+h)-F(x) / h to fdash=F(x+h)-F(x-h) / 2h
This is a good compromise between tersity of code and precision.

Also in you notes what is the symbol in 1st derivatives

fdash = m1 + unknown_symbol (dx^2)The symbol dx is the step-size, often called h.

Also how do you access g++ documentation/help files?You can look at the manual pages of g++ by typing the command 'man g++'. There is also a good documentation available online from GNU. I'll get the link tomorrow. I also have this book (but it does not contain any more than the above mentioned docs): GCC book (http://www.amazon.com/Definitive-Guide-GCC-Second/dp/1590595858/ref=pd_bbs_sr_1?ie=UTF8&s=books&qid=1217189173&sr=1-1)

P.S do you think it is a good Idea to change the title of this thread to "Numerical programing" or something like that?We shouold keep going in this thread for this assignment. It is not unusual to have a thread go for 50 or more entries over the course of a few weeks. If we had done anything wrong, the mods would have commented. After a while, certain threads like this one don't draw the full attention of the users. So if you have an entirely different question, then you might want to post an additional thread. Also I might mention that I am not often on the board because I have very many other topics of work. But I will continue checking your progress on this thread.

Feel free to post additional questions as you progress.

Sincerely, Chris.

[Nighteyes]

Just posting a little update.

Read another 400 pages of my c++ book. read:

8) operator overloading
9)Inheritance
10)Pointers
11)Virtual Functions
12) Streams and Files
13) Multifile Programs

I now understand why you made the namespace's.

So now I made the Ashley program so it imported into the Newton_R method solver and finds the zero.

At the moment there is a big error in the ashley program as can be seen in the pic "ah.jpg" this is not due to the code though this is due to the variables A, B, C and D which I will renter tomorrow. Otherwise I am quite pleased with my first attempt at a multifile program. Though I think I am going to do alot of changes to how the program will work in the end.

I am thinking have a input file which reads in different omega values and initial guesses and then outputs the final k value along with the inputed omega value.

On a side point to make the Newton_R Method program an independent module which can be used to solve all sorts of equations would it be a good Idea to make it in to a class? or at least a function that takes an input of the function that needs to be solved and an initial guess....

Thanks Craig

[dude_1967]

Hi Craig,

Sorry about the late response. I was unavailable for a while.

You should probably try to have a separate cpp file for the Ashley stuff. By putting the static functions right in the header file, you have essentially just included one cpp file into another. You need to make a separate header file for the Ashley stuff which includes only the function declarations in the namespace---and only those needed by other modules.

I am having difficulty understanding you implementation of Newton-Raphson. I am unable to see where you re-evaluate the function and its derivative using the new midpoint. I am unsure, but I beleive that something might be wrong here. But I have written a sample included below which might help you clear this up.

As to your question about a generic implementation of Newton-Raphson, I would definitely consider this to be the right thing to do. However this is often simply implemented as a procedure, not a class. You have an additional complication with your independent function GA::Eigen. It takes two input parameters, the second one named omega. Most simple, generic implementations of Newton-Raphson procedures deal with functions with only one independent variable, say x.

I have written down an implementation of Newton-Raphson using a third order derivative formula in the test program shown below. It is based on the algorithm shown in the book "Numerical Recipes in C++". You should look carefully how a function pointer is passed as an input parameter to the Newton-Raphson method as well as to the derivative algorithm. This technique is essential if you want to create any kind of generic algorithm for these kinds of things.

We could modify this technique to accomidate your function GA::Eigen(x, omega) if you would like. You could also try to do this on your own and get back to me with any problems.

Good luck with your further progress.

Sincerely, Chris.