pow: DOMAIN error

[ellisrn]

Hi,

I'm trying to run this short code and when it prints out i get about 20 of these "pow: DOMAIN error" messages, and then my results. Can anyone shed any light on this. I'm not acutally a programer, but got thrust into this in a project for another class, so I'm extremely clueless. I think it might have something to do with the output file. I don't really know what that is, and I don't even know if what I have in there is legitimate.

Thanks!

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <fstream.h>

//Global Variables

const int nvar   =3;
const int nmom   =3;

double delta;


double A[nvar+1][nvar+1];
double B[nvar+1][nvar+1];
double C[nvar+1][nvar+1];
double beni[nvar+1];
double benij[nvar+1][nvar+1];

double v[nvar+1][nmom+1];     //in delme norberg
double dydx[nvar+1][nmom+1];

double mu[nvar+1][nvar+1];

//Program Headers

void derivs(double t, double yt[nvar+1][nmom+1], double dyt[nvar+1][nmom+1]);
double fact (int n);
double comb (int i,int j);
double power (double i,double j);


//Specify Output File

ofstream outfile

("h:/bcppexe/derivs/output/test.txt", ios::out);   //how do you direct an output file?




main()
{
    double delta = 0.05;





   A[1][2] = 0.0004;
   A[1][3] = 0.0005;
   A[2][3] = 0.0005;
   A[2][1] = 0;
   A[3][1] = 0;
   A[3][2] = 0;

   B[1][2] = 0.0000034674;
   B[1][3] = 0.000075858;
   B[2][3] = 0.000075858;
   B[2][1] = 0;
   B[3][1] = 0;
   B[3][2] = 0;

   C[1][2] = 0.06;
   C[1][3] = 0.038;
   C[2][3] = 0.038;
   C[2][1] = 0;
   C[3][1] = 0;
   C[3][2] = 0;

   A[1][1] = 0;
   B[1][1] = 0;
   C[1][1] = 0;
   A[2][2] = 0;
   B[2][2] = 0;
   C[2][2] = 0;
   A[3][3] = 0;
   B[3][3] = 0;
   C[3][3] = 0;

    beni[1] = 0;
   beni[2] = 0;
   beni[3] = 0;

   benij[1][3] = 1;
   benij[2][3] = 1;
   benij[1][2] = 0;
   benij[2][1] = 0;
   benij[3][1] = 0;
   benij[3][2] = 0;
   benij[1][1] = 0;
   benij[2][2] = 0;
   benij[3][3] = 0;


   for(int i=0; i <= nvar; i++)
   {
    for(int j=0; j <= nmom; j++)
      {
        v[i][j] = 0;         //Initializing the matrix of yt values to zero, and setting the boundary condition
         dydx[i][j] = 0;      //Initializing the matrix of the ODE to zero
         v[i][0] = 1;         //Setting the boundary condition for the 0th moment
      }
   }


   derivs(60, v, dydx);








  return 0;


}



void derivs(double t, double yt[nvar+1][nmom+1], double dyt[nvar+1][nmom+1])
{

   int i;
   int j;
   int r;
   int q;


   double summation[nvar +1][nmom +1];
   double mu_start[nvar+1];
   double sum_part[nvar+1][nmom+1];



    double mu[nvar+1][nvar+1];



        for(i=1; i <= nvar; i++)
    {
        for(j=1; j <= nvar; j++)
        {
                mu[i][j] = A[i][j] + B[i][j] * pow(10, (C[i][j]*t));
        }
    }





   for(i=1; i <= nvar; i++)
   {
      mu_start[i] = 0;

    for(j=1; j <= nvar; j++)
      {
         mu_start[i] += mu[i][j];
      }
   }




  for(i=1; i <= nvar; i++)         //sums up the second part of the summation in Norberg's equation
  {

      for(q=1; q <= nmom; q++)
      {

            sum_part[i][q] = 0;

            for(j=1; j <= nvar; j++)
        {

            for(r=0; r <= q; r++)
            {

               sum_part[i][q] += comb(q,r) * pow(benij[i][j], r) * yt[j][q-r];
            }
         }
      }
  }
















   for(i=1; i <= nvar; i++)             //multiplies and stores the total summation in Norberg's equation
   {

      for(q=1; q <= nmom; q++)
      {

        summation[i][q] = mu_start[i] * sum_part[i][q];
      }
   }







    for(i=1; i <= nvar; i++)
   {
    for(q=1; q <= nmom; q++)
      {
         yt[i][q] = 0;

         yt[i][0] = 1;



            dyt[i][q] = ((((q * delta) + mu_start[i]) * yt[i][q]) - (q * beni[i] * yt[i][q-1]) - (summation[i][q]));
         }
      }








   for(i=1; i <= nvar; i++)
   {
    for(q=1; q <= nmom; q++)
      {
        cout << dyt[i][q] << endl;
      }
   }


   char w;
   cin>>w;














}














double fact (int n)
{
        if (n < 0) return 0;
    double f = 1;
    while (n > 1) f *= n--;

    return f;
}




double comb (int i,int j)
{
    double c = 1;
        if (j==0) return c;
    else      c = (fact (i)) / (fact (j) * fact (i-j));

    return c;
}





 double power (double i,double j)
{
    double p=1;
    if (j == 0) return 1;
    else if (i==0) return 0;
    else p = pow(i,j);

    return p;
}

[Dave Sinkula]

>Can anyone shed any light on this.

The pow functions compute x raised to the power y. A domain error occurs if x is finite and negative and y is finite and not an integer value. A domain error may occur if x is zero and y is less than or equal to zero. A range error may occur.

for ( int i=0; i <= nvar; i++ )
{
	 for ( int j=0; j <= nmom; j++ )

You should really learn how to index arrays properly. The usual idiom is as follows.

double array[3][4];
for(int i = 0; i < 3; ++i)
{
   for(int j = 0; j < 4; ++j)
   {
	 /* ... */
   }
}

That is, an array of size N may be indexed from 0 to N-1.