What is causing the incorrect rendering of the Feigenbaum tree in C with OpenGL?

  • Thread starter computerex
  • Start date
  • Tags
    Tree
In summary, the conversation discusses an attempt to render the bifurcation diagram of the logistic difference equation using C with OpenGL. The program produces output, but the structure does not resemble the desired result. The conversation suggests checking the values being plotted and plotting all different values in order to create the complete graphic. Some implementations plot (g(x), x) where g(x) = X(n+1) = k*Xn*(Xn-1), while others simply have the k being tested as the x value. The conversation also notes a potential issue with the definition of NUM_OF_ITT.
  • #1
computerex
68
0
Hello. I am attempting to render the bifurcation diagram of the logistic difference equation in C with OpenGL. I think the math behind it is simple enough, but for some reason I am ending up with a structure that resembles the tree but doesn't look quite like it. I would really appreciate it if someone could point out where I am going wrong. Here is how it should look like:

http://brain.cc.kogakuin.ac.jp/~kanamaru/Chaos/e/BifArea/
http://planetmath.org/encyclopedia/FeigenbaumFractal.html

The following program produces the output attached to this post. As you can see, you can distinguish that the structure seems to be correct.

Code: 
#include <stdio.h>
#include <stdlib.h>
#include <gl/glut.h>
#include <windows.h>
#include <math.h>

#define NUM_OF_ITT  100.0 // number of itterations
#define NUM_OF_SETS 1000  //  number of ks that will be tested

double xstart = 0.1;      // x0 
double dataPoints[NUM_OF_SETS][(int)NUM_OF_ITT][2]; //

void calculatePlot(double k, unsigned int inx)
{
	double result = k*xstart*(1.0-xstart);
	dataPoints[inx][0][0]=k;
	dataPoints[inx][0][1]=result;
	for ( unsigned int i = 1; i < NUM_OF_ITT; i++ )
	{
		result = k * result*(1.0-result);
		dataPoints[inx][i][0]=k;
		dataPoints[inx][i][1]=result;
	}
}

void handleKeyPress(unsigned char key, int x, int y)
{
	switch (key)
	{
	case 27:
		exit(0);	
	}
}
void initRendering()
{
	glEnable(GL_DEPTH_TEST);
	glLineWidth(1.0);
	double k = 1;
	double kstep = (4.0-k)/(double)NUM_OF_SETS;
	for ( unsigned int i = 0; i < NUM_OF_SETS; i++ )
	{
		calculatePlot(k, i);
		k+=kstep;
	}
}

void handleResize(int q, int z)
{
	glViewport(0, 0, q, z);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(45.0, (double)q / (double)z, 1.0, 200.0);
}

void drawScene()
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);

	glLoadIdentity();
	glTranslatef(-3, 0, -5);
	glColor3f(0.0, 1.0, 0.0);

	// not rendering the transient states
	glBegin(GL_LINE_STRIP);
	for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
	{
		glVertex2f(dataPoints[set][0][0], dataPoints[set][0][1]);
	}
	glEnd();
	glBegin(GL_LINE_STRIP);
	for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
	{
		glVertex2f(dataPoints[set][(int)NUM_OF_ITT-1][0], dataPoints[set][(int)NUM_OF_ITT-1][1]);
	}
	glEnd();
	glutSwapBuffers();
}

void update(int value)
{
	glutTimerFunc(10, update, 0);
}

void idlefunc()
{
	glutPostRedisplay();
}
void mouse(int button, int state, int x, int y)
{
}
int main(int argc, char **argv)
{
	glutInit(&argc, argv);	
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB |GLUT_DEPTH);
	glutInitWindowSize(400, 400);
	glutCreateWindow("Logistcs");
	initRendering();
	glutDisplayFunc(drawScene);
	glutIdleFunc(idlefunc);
	glutKeyboardFunc(handleKeyPress);
	glutReshapeFunc(handleResize);
	glutMouseFunc(mouse);
	glutTimerFunc(100, update, 0);
	glutMainLoop();
	return 0;
}
 

Attachments

  • feigen.png
    feigen.png
    609 bytes · Views: 520
Last edited by a moderator:
Technology news on Phys.org
  • #2
I think you are only displaying the data point for the last iteration.

Print out say the last few iteration values in the datapoint array, in each call to calculateplot(). When k is below 3, you will find they are all the same. Then, they start to flip between 2 diffeent values. Then they go round a cycle of 4 different values, then 8, etc. Later on, they are chaotic. Later still they will repeat in cycles with other lengths, like 6 and 3.

You need to plot all the different values to make the complete graphic. One way is to calculate say 1000 iterations, to wash out the effect of your choice of a starting value, and then plot say the next 100 iterations. Many of the 100 plotted points will lie on top of each other, but that doesn't matter.
 
  • #3
AlephZero said:
I think you are only displaying the data point for the last iteration.

Hello. Thank you for the reply.
Take a look at this:

Code: 
double dataPoints[NUM_OF_SETS][(int)NUM_OF_ITT][2];

This multidimensional array contains the x-y coordinate data for each k that is going to be tested. The number of different k's tested will be NUM_OF_SETS k's. The difference of some k and the next k tested will be 3/NUM_OF_ITT. Here's the routine that precomputes the data for better rendering times. I think the problem may lie here:

Code: 
void calculatePlot(double k, unsigned int inx)
{
	double result = k*xstart*(1.0-xstart);
	dataPoints[inx][0][0]=k;
	dataPoints[inx][0][1]=result;
	for ( unsigned int i = 1; i < NUM_OF_ITT; i++ )
	{
		result = k * result*(1.0-result);
		dataPoints[inx][i][0]=k;
		dataPoints[inx][i][1]=result;
	}
}

0 <= inx < NUM_OF_SETS. Note how the x coordinate is always set to the k being tested. I am not sure whether this is the right way of doing it. Some implementations plot (g(x), x) where g(x) = X(n+1) = k*Xn*(Xn-1). Others simply have the k being tested as the x value. Other still use (x, g(x)).

This is the rendering procedure, which is relatively simple. Note, right now for each k being tested the x value for each point is k itself, and only y varies. And so you would expect to see a vertical line for each k, with the points between dataPoints[set][0] and dataPoints[set][NUM_OF_ITT-1] being transient states. Therefore those points are not rendered here. However, that data is still available and can be rendered. Indeed originally it was rendered, but I saw only vertical columns connecting the top outline and the bottom outline. So I simply decided to omit those points from rendering.

The mechanism of rendering is simple: We loop through all the sets (the data for all the k's tested) and render a line strip using the first point calculated. Then the same procedure is repeated, rendering the last data point in the array for each set.

Code: 
// not rendering the transient states
glBegin(GL_LINE_STRIP);
for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
{
	glVertex2f(dataPoints[set][0][0], dataPoints[set][0][1]);
}
glEnd();

glBegin(GL_LINE_STRIP);
for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
{
	glVertex2f(dataPoints[set][(int)NUM_OF_ITT-1][0], dataPoints[set][(int)NUM_OF_ITT-][1]);
}
glEnd();

I don't know what I am doing wrong. I am sorry, but I really have no background education of dynamic systems (being a high school senior). I just love math/realtime computer graphics. Further help would be much appreciated.
 
  • #4
Minor points
Why is NUM_OF_ITT defined as a floating point value?
Code: 
#define NUM_OF_ITT  100.0 // number of itterations
This doesn't make sense, especially since you have to cast it to int when you declare your dataPoints array.

BTW, iteration doesn't have any double t's in it, which makes NUM_OF_ITT a bit more difficult to understand than it needs to be.
 
  • #5
Mark44 said:
Minor points
Why is NUM_OF_ITT defined as a floating point value?
Code: 
#define NUM_OF_ITT  100.0 // number of itterations
This doesn't make sense, especially since you have to cast it to int when you declare your dataPoints array.

BTW, iteration doesn't have any double t's in it, which makes NUM_OF_ITT a bit more difficult to understand than it needs to be.

I thought I'd have to do some arithmetic later that'd involve casting the constant to a floating point later on, and I was always bad at spelling ;)

In any case, you can be rest assured that neither of the aforementioned minor programming inconveniences are the cause of the issue. Do you have any other ideas :) ?
 
  • #6
I wasn't able to open the first link you posted. I was able to open the planetmath.org link, however.

Another (and clearer, IMO) explanation that seems closer to what you're trying to do is here: http://mathforum.org/advanced/robertd/bifurcation.html

The bifurcation graph represents the results of repeated calculating [itex]\mu[/itex]x(1 - x) for values of [itex]\mu[/itex] and x. (You're using k for [itex]\mu[/itex] and y for x.)

It works by picking a value of [itex]\mu[/itex], and a starting value of x. The next value of x is obtained by substituting it into the expression [itex]\mu[/itex]x(1 - x). For some values of [itex]\mu[/itex] the iterated values converge to a specific number.

For example, if [itex]\mu[/itex] is 2 and we start with x = .25 and calculate 2*.25*.75, the next x value will be .375. Putting this value in for x, the next x value will be .46875. After a few steps you get values that keep repeating, .5, .5, .5. So one point on the graph you're trying to get will be (2, .5). Similarly, when [itex]\mu[/itex] is 2.5, the iterations settle down to .6, so another point on the graph is (2.5, .6).

For values of [itex]\mu[/itex] from 1 to about 3, iteration produces a single value. For larger values of [itex]\mu[/itex], iteration can produce two values, then four values, then eight values, and so on, and that's what gives you the bifurcations, or forks in the diagram. When [itex]\mu[/itex] is around 3.5, you start to get cycles of four different iterated values, so the bifurcation graph will show four different points for a while.

Now maybe you know this or maybe not. Your code doesn't contain may comments, so it's hard to tell what you're doing, particularly with your three dimensional array. It would be nice to have a comment saying what each index in that array represents.
 
  • #7
I am sorry :( I know it can be a nightmare trying to interpreted undocumented code, so I will try to explain to the best of my abilities what I am doing in the code:

Code: 
#define NUM_OF_IT  100       // number of iterations
#define NUM_OF_SETS 1000  //  number of ks that will be tested

double xstart = 0.1;      // x0 
double dataPoints[NUM_OF_SETS][(int)NUM_OF_IT][2]; //

The basic algorithm I am following is:

1) Pick some starting x value that is below 1.0. This is xstart in the program (0.1).
2) Pick some K <= 1 <= 4
3) Iterate using a for loop NUM_OF_IT times. In each iteration, use the logistics difference equation to calculate the y value, and store it. This y value will be used in the next iteration. For the first iteration, the y value is calculated as: y = K*xstart*(1-xstart). This is stored in the dataPoints array, which has this structure:

dataPoints[index for data for some unique K][coordinate pair calculated in the ith iteration] = [k value (constant for this index), y value]

4) Render to the screen

This will become clearer when taking a look at the following code:

Code: 
// This is where the iteration takes place
// inx is the index for the dataPoints array where the data for this unique k will be stored
void calculatePlot(double k, unsigned int inx)
{
        // this is the base case
        // all subsequent values of y will be calculated using this
	double result = k*xstart*(1.0-xstart);

        // dataPoints[inx][i][0] will [i]always[/i] contain k
	dataPoints[inx][0][0]=k;
	dataPoints[inx][0][1]=result;
	
        for ( unsigned int i = 1; i < NUM_OF_IT; i++ )
	{
                // the previous result is used to calculate the new value for y
		result = k * result*(1.0-result);
		dataPoints[inx][i][0]=k;
		dataPoints[inx][i][1]=result; // y is stored in the 1st element of the 3rd 
                                                      //  dimension of dataPoints array
	}
}

Code: 
void initRendering()
{
        // standard OpenGl initialization. 
	glEnable(GL_DEPTH_TEST);
	glLineWidth(1.0);

        // Here is where the dataPoints array is actually populated with data. 
        // the first k we use is 1. Then we increase k by a constant quantity 
        // so the final k is <= 4

	double k = 1;
	double kstep = (4.0-k)/(double)NUM_OF_SETS;
	for ( unsigned int i = 0; i < NUM_OF_SETS; i++ )
	{
                // The data for this particular k is stored in the ith element of the 
                // first dimension of the dataPoints array
		calculatePlot(k, i);
                // advance k by the constant quantity
		k+=kstep;
	}
}

Code: 
        // From what I have seen, some implementations only render the 
        // x-y values of the first iteration and the final iteration. This is what I have
        // done here. 
	glBegin(GL_LINE_STRIP);
        // pseudocode: "for each value of K, render the point calculated in the first iteration"
	for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
	{
		glVertex2f(dataPoints[set][0][0], dataPoints[set][0][1]);
	}
	glEnd();

        // pseudocode: "for each value of K, render the point calculated in the last iteration"
	glBegin(GL_LINE_STRIP);
	for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
	{
		glVertex2f(dataPoints[set][(int)NUM_OF_IT-1][0], dataPoints[set
[(int)NUM_OF_IT-1][1]);
	}
	glEnd();

I hope this was enlightening. I am experimenting locally, hopefully I'll come up with a working solution soon. In the meantime, and help will be very much appreciated.
 
  • #8
I am the biggest idiot in the world. As opposed to drawing lines, I should have simply drawn points. The structure seemed correct because it was. There was simply extra information there. I hope I am making sense. Thanks for the help all. Attached picture + the working source below.

Code: 
#include <stdio.h>
#include <stdlib.h>
#include <gl/glut.h>
#include <windows.h>
#include <math.h>

#define NUM_OF_ITT  1000.0 // number of itterations
#define NUM_OF_SETS 1000  //  number of ks that will be tested

double xstart = 0.1;      // x0 
double dataPoints[NUM_OF_SETS][(int)NUM_OF_ITT][2]; //

void calculatePlot(double k, unsigned int inx)
{
	double result = k*xstart*(1.0-xstart);
	dataPoints[inx][0][0]=k;
	dataPoints[inx][0][1]=result;
	for ( unsigned int i = 1; i < NUM_OF_ITT; i++ )
	{
		result = k * result*(1.0-result);
		dataPoints[inx][i][0]=k;
		dataPoints[inx][i][1]=result;
	}
}

void handleKeyPress(unsigned char key, int x, int y)
{
	switch (key)
	{
	case 27:
		exit(0);	
	}
}
void initRendering()
{
	glEnable(GL_DEPTH_TEST);
	glLineWidth(1.0);
	glPointSize(1.0);
	double k = 1;
	double kstep = (4.0-k)/(double)NUM_OF_SETS;
	for ( unsigned int i = 0; i < NUM_OF_SETS; i++ )
	{
		calculatePlot(k, i);
		k+=kstep;
	}
}

void handleResize(int q, int z)
{
	glViewport(0, 0, q, z);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(45.0, (double)q / (double)z, 1.0, 200.0);
}

void drawScene()
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);

	glLoadIdentity();
	glTranslatef(-3, 0, -5);
	glColor3f(0.0, 1.0, 0.0);

	for ( unsigned int set = 0; set < NUM_OF_SETS; set++ )
	{
		glBegin(GL_POINTS);
		for ( unsigned int i = 0; i < NUM_OF_ITT; i++ )
		{
			glVertex2f(dataPoints[set][i][0], dataPoints[set][i][1]);
		}
		glEnd();
	}
	glutSwapBuffers();
}

void update(int value)
{
	glutTimerFunc(10, update, 0);
}

void idlefunc()
{
	glutPostRedisplay();
}
void mouse(int button, int state, int x, int y)
{
}
int main(int argc, char **argv)
{
	glutInit(&argc, argv);	
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB |GLUT_DEPTH);
	glutInitWindowSize(400, 400);
	glutCreateWindow("Logistcs");
	initRendering();
	glutDisplayFunc(drawScene);
	glutIdleFunc(idlefunc);
	glutKeyboardFunc(handleKeyPress);
	glutReshapeFunc(handleResize);
	glutMouseFunc(mouse);
	glutTimerFunc(100, update, 0);
	glutMainLoop();
	return 0;
}
 

Attachments

  • bifurcation.png
    bifurcation.png
    2.5 KB · Views: 521
  • #9
Well, it's not quite right yet, because you have those strange lines from the bottom left hand corner.

They are the points when you "start up" the iteration from x = 0.1, until it settles down.

If you change 0.1 to a differrent value (say 0.9) you will get a different pattern of those lines.

To get rid of them completely, change the line in drawscene() from

for ( unsigned int i = 0; i < NUM_OF_ITT; i++ )

to something like

for ( unsigned int i = NUM_OF_ITT/2; i < NUM_OF_ITT; i++ )

to skip over the first points you calculate.
 
  • #10
Mark44 said:
I wasn't able to open the first link you posted.
I thnk Japan is off-line right now :rolleyes:
 
  • #11
AlephZero said:
I thnk Japan is off-line right now :rolleyes:
I didn't look at the link that closely to notice that it was a site in Japan. They have other problems right now...
 
  • #12
The strange thing is, the site was working perfectly when I made the first post. It stopped responding very soon after.
 

1. What is the Feigenbaum tree?

The Feigenbaum tree is a mathematical object that shows the universality of the period-doubling route to chaos in a simple, one-dimensional dynamical system called the logistic map. It was discovered by mathematician Mitchell Feigenbaum in the 1970s.

2. How is the Feigenbaum tree rendered?

The Feigenbaum tree is rendered using a process called iteration. This involves repeatedly applying a mathematical formula to each point on a number line, and then plotting the resulting points on a graph. The process is repeated with different starting values, and the points are connected to create the Feigenbaum tree.

3. What does the Feigenbaum tree look like?

The Feigenbaum tree is a fractal object that has a self-similar structure at different scales. It starts with a main trunk that branches out into smaller and smaller branches, creating a tree-like shape. The branches are not evenly spaced, but instead follow a specific mathematical pattern.

4. What is the significance of the Feigenbaum tree?

The Feigenbaum tree is significant because it provides insight into the behavior of chaotic systems. It shows that even in a seemingly simple and predictable system, small changes in initial conditions can lead to unpredictable and chaotic behavior. This has important implications in fields such as physics, biology, and economics.

5. How is the Feigenbaum tree related to the concept of universality?

The Feigenbaum tree is an example of universality, which is the idea that certain patterns or behaviors are repeated across different systems. In this case, the period-doubling route to chaos is a universal phenomenon that can be observed in many different systems, from the growth of bacterial colonies to the behavior of the stock market.

Similar threads

Trying to compile a source code in C
    • Programming and Computer Science
Replies
5
Views
1K
Need help in a Basket Ball game(Projectile Motion / motion under Gravitation)
    • Programming and Computer Science
Replies
1
Views
7K
Fixing OpenGL App: Black Background & White Triangle
    • Computing and Technology
Replies
2
Views
2K

Hot Threads

Recent Insights

Back
Top