Hey there^^
Just a quick sort out hopefully, ive been hacking at this for a good while and cant figure out whats wrong for thie life of me! I have a program that uses calculates the maximum deflection of a beam under load (with some assumptions were allowed); and im now trying to add to that program a for loop that repeats it 1000 times (simple enough), and a function that calculates random stiffness values for each repetition. Trouble is, the function is giving me this funny error message :/
Anyhoo, here it is; the errors are generating on lines 27 and 31.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define L 0.8
#define F 10000
#define E 1.36E+10 //defined constants.
#define s_dev E+9
#define delta_x 0.016
double RandInvNormal(double , double); //random stiffness function header.
int main(int argc, char *argv[])
{
double I[51], C[51], y[51], z[50]; //double precision arrays.
float percentage_error; //counter for error convergence.
int i, x; //increment counter for array assignment.
union { unsigned int i[2]; time_t d;} t;
t.d=time(NULL);
srand(t.i[0]);
for(x = 1; x < 1000; x++) //repeats entire algorithm for 1000 iterations.
{
for (i = 0; i < 51; i++) //for loop assigns array values.
{
I[i] = 6.833E-07 *exp ( -19.02 *pow( ((delta_x*i) - (L/2)), 2 )); //calculates and assigns I values into array.
if((delta_x*i)<= L/2) //calculates and assigns 1/R values into C[]array, using I values.
{
C[i]=F*(delta_x*i)/(2.0*RandInvNormal(&E, &s_dev)*I[i]);
}
else
{
C[i]=F*(L-(delta_x*i))/(2.0*RandInvNormal(&E, &s_dev)*I[i]);
}
y[i] = 0; //sets initial guess for deflection to zero.
}
percentage_error = 1; //resets error inside greater for loop.
while (percentage_error >= 0.001) //while loop for testing change in max deflection.
//continues if change is greater than 0.001%.
{
for (i = 1; i < 50; i++) //loop iterates previous values.
{
z[i] = (y[i+1] + y[i-1] + (C[i] *pow(delta_x, 2)))/2;
//simply supported ends - zero vertical deflection.
z[0]=0;
}
percentage_error = (z[25] - y[25]) *(100/z[25]); //test while loop condition.
for (i = 1; i < 50; i++) //loop appends y-array with new set of iterations for iterative loop.
{
y[i] = z[i];
}
}
printf("\nY(max) = %lf\n", y[25]); //displays maximum deflection of beam *assumed to be at mid-point*
}
}
//------------------------------FUNCTIONS---------------------------------
double RandInvNormal(double mean, double stdev)If you'd like to run the program and test it out, the function code is:
{
// Define constants to curve fit for integral
const double A1 = -3.969683028665376e+01;
const double A2 = 2.209460984245205e+02;
const double A3 = -2.759285104469687e+02;
const double A4 = 1.383577518672690e+02;
const double A5 = -3.066479806614716e+01;
const double A6 = 2.506628277459239e+00;
const double B1 = -5.447609879822406e+01;
const double B2 = 1.615858368580409e+02;
const double B3 = -1.556989798598866e+02;
const double B4 = 6.680131188771972e+01;
const double B5 = -1.328068155288572e+01;
const double C1 = -7.784894002430293e-03;
const double C2 = -3.223964580411365e-01;
const double C3 = -2.400758277161838e+00;
const double C4 = -2.549732539343734e+00;
const double C5 = 4.374664141464968e+00;
const double C6 = 2.938163982698783e+00;
const double D1 = 7.784695709041462e-03;
const double D2 = 3.224671290700398e-01;
const double D3 = 2.445134137142996e+00;
const double D4 = 3.754408661907416e+00;
const double P_LOW = 0.02425;
const double P_HIGH = 0.97575; // P_high = 1 - p_low
double p, q, r, x;
do
{
p = (float)(rand())/RAND_MAX; // Get a random number between 0 and 1
// rand() returns an int between 0 and RAND_MAX
} while (p <= 0 || p >= 1); // Repeat making sure p is not 0 or 1
if (p > 0 && p < P_LOW)
{
q = sqrt(-2*log(p));
x = (((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6) / ((((D1*q+D2)*q+D3)*q+D4)*q+1);
} else if (p >= P_LOW && p <= P_HIGH)
{
q = p - 0.5;
r = q*q;
x = (((((A1*r+A2)*r+A3)*r+A4)*r+A5)*r+A6)*q /(((((B1*r+B2)*r+B3)*r+B4)*r+B5)*r+1);
} else {
q = sqrt(-2*log(1-p));
x = -(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6) / ((((D1*q+D2)*q+D3)*q+D4)*q+1);
}
return (mean + x*stdev);
}Thanks for any help :)
