[sabago]

Dear OpenFoamers!

I'm struggling with getting my root finding (bisection) method to work although the pseudo-code gives me the correct root in an ordinary C++ compiler.

Below is what I have in OF

{
forAll(x1, celli)
forAll(x2, celli)
forAll(x3, celli)
forAll(f1, celli)
forAll(f3, celli)
forAll(root, celli)
{
 do
    {
       if(count == iter)
       {
           break;
       }
 
        x3[celli] = (x1[celli] + x2[celli])/2;
Info<<x3<<endl;

	f1 = (Foam::exp(8.1*x1))-(Foam::exp(-2.7*x1))-(j/i);
Info<<f1<<endl;

	f3 = (Foam::exp(8.1*x3))-(Foam::exp(-2.7*x3))-(j/i);
Info<<f3<<endl;

	 if((f1[celli] *f3[celli])  < scalar(0.0) )
        {
            x2[celli] = x3[celli];
        }
        else
        {
            x1[celli] = x3[celli];
        }
        count =count+1;
    }
    while ((mag(x1[celli] - x2[celli]) < scalar(1e-12))||(f3[celli]==0));

  root[celli] =x3[celli];
Info<<root<<endl;
}
}

where j=-2e6 and i is calculated elsewhere in OF and is about 5e9.
my x1 is set to 0 and x2 to -1 as the root is about -3e-5; The first step works fine in OpenFOAM with iter=100, I can see in the terminal that x3 is indeed -0.5 but I get a root of 0! Is it because the root is so small?

Below is the pseudo-code that gives me a root of -3.3e-5 in C++ (if needed for comparison)

#include <iostream>
//#include<conio.h> DOS hearder, use curses.h instead
#include <ncurses.h>
#include<math.h>
#include <stdlib.h>
#include <math.h>  
//#include <chplot.h>
using namespace std;

 
int main()
{

	float x1 = 0;
	float x2 =-1;
	float x3;
	int count =0;
	int iter = 20;
    do
    {
       if(count == iter)
       {
           break;
       }
 
        x3 = (x1 + x2)/2;
 
       cout <<"x1=" << x1 <<" | x2="<< x2 <<" | x3=" << x3 << endl << endl;
 
        //float temp1 = f(x1);
        //float temp2 = f(x3);
        float f1 =exp(8.1*x1)-exp(-2.7*x1) + 0.00037;
        float f3 =exp(8.1*x3)-exp(-2.7*x3) + 0.00037;
         //cout <<"f1=" << f1 <<endl;
        //   cout <<"f3=" << f3 <<endl;
        if((f1) * (f3) < 0 )
        {
            x2 = x3;
        }
        else
        {
            x1 = x3;
        }
        count++;
    }
    while ((abs(x1 - x2) < 0.0000001) || (exp(8.1*x3)-exp(-2.7*x3) + 0.00037 == 0 ));
 cout << "The root value is " << x3 << endl;
}

Any help will be much appreciated.

Best,
Sandra

[sabago]

Dear all,
I managed to fix the code (please see below). It turns out that I didn't need the multiple forAll caluses (adding them all in significantly slows down the program) and you need at least one forAll declaration otherwise the if statements give errors. I have a question, however, can I have a "while" loop in OF? I think the issue was that OF has a different style of "while" loop. Does anyone have an idea? I checked this link: http://www.openfoam.org/contrib/code-style.php to see what a "while" loop should look like but there doesn't seem to be anything on while loops.

Anyway, below is the revised code (iter=100):

forAll(x3, celli)
{
for (count=0; count<iter; count++) 
{
        x3 = (x1 + x2)/2;

	f1 = (Foam::exp(8.1*x1))-(Foam::exp(-2.7*x1))-(j/i);

	f3 = (Foam::exp(8.1*x3))-(Foam::exp(-2.7*x3))-(j/i);
	
	 if(neg(f1[celli] *f3[celli]))
        {
            x2 = x3;
        }
        else if((mag(x1[celli] - x2[celli]) < 1e-7)||(f3[celli]==0))
        {
            root =x3;
	    break;
        }
	else
	{
	x1= x3;
	}
}
     
Info<<"Root is"<<endl;
Info<<root<<endl;
}

Hope this helps somebody else.

Best,
Sandra

Originally Posted by sabago

Dear OpenFoamers!

I'm struggling with getting my root finding (bisection) method to work although the pseudo-code gives me the correct root in an ordinary C++ compiler.

Below is what I have in OF

{
forAll(x1, celli)
forAll(x2, celli)
forAll(x3, celli)
forAll(f1, celli)
forAll(f3, celli)
forAll(root, celli)
{
 do
    {
       if(count == iter)
       {
           break;
       }
 
        x3[celli] = (x1[celli] + x2[celli])/2;
Info<<x3<<endl;

	f1 = (Foam::exp(8.1*x1))-(Foam::exp(-2.7*x1))-(j/i);
Info<<f1<<endl;

	f3 = (Foam::exp(8.1*x3))-(Foam::exp(-2.7*x3))-(j/i);
Info<<f3<<endl;

	 if((f1[celli] *f3[celli])  < scalar(0.0) )
        {
            x2[celli] = x3[celli];
        }
        else
        {
            x1[celli] = x3[celli];
        }
        count =count+1;
    }
    while ((mag(x1[celli] - x2[celli]) < scalar(1e-12))||(f3[celli]==0));

  root[celli] =x3[celli];
Info<<root<<endl;
}
}

where j=-2e6 and i is calculated elsewhere in OF and is about 5e9.
my x1 is set to 0 and x2 to -1 as the root is about -3e-5; The first step works fine in OpenFOAM with iter=100, I can see in the terminal that x3 is indeed -0.5 but I get a root of 0! Is it because the root is so small?

Below is the pseudo-code that gives me a root of -3.3e-5 in C++ (if needed for comparison)

#include <iostream>
//#include<conio.h> DOS hearder, use curses.h instead
#include <ncurses.h>
#include<math.h>
#include <stdlib.h>
#include <math.h>  
//#include <chplot.h>
using namespace std;

 
int main()
{

	float x1 = 0;
	float x2 =-1;
	float x3;
	int count =0;
	int iter = 20;
    do
    {
       if(count == iter)
       {
           break;
       }
 
        x3 = (x1 + x2)/2;
 
       cout <<"x1=" << x1 <<" | x2="<< x2 <<" | x3=" << x3 << endl << endl;
 
        //float temp1 = f(x1);
        //float temp2 = f(x3);
        float f1 =exp(8.1*x1)-exp(-2.7*x1) + 0.00037;
        float f3 =exp(8.1*x3)-exp(-2.7*x3) + 0.00037;
         //cout <<"f1=" << f1 <<endl;
        //   cout <<"f3=" << f3 <<endl;
        if((f1) * (f3) < 0 )
        {
            x2 = x3;
        }
        else
        {
            x1 = x3;
        }
        count++;
    }
    while ((abs(x1 - x2) < 0.0000001) || (exp(8.1*x3)-exp(-2.7*x3) + 0.00037 == 0 ));
 cout << "The root value is " << x3 << endl;
}

Any help will be much appreciated.

Best,
Sandra