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;
}