| T O P I C R E V I E W |
| espenhjo |
Posted - 10/15/2003 : 08:11:26 AM
Dear reader,
I have a problem consisting of an NLSF of raw data to a function that has transcendental equations. My data can be fitted to the equation:
I guess that the equation converges quite fast, so that n to 100 would be sufficient. The xis and zetas are roots to the transcendental equations
I want to fit the data in Origin 7 using the NLSF tool with Origin C.
The definition in NLSF would look something like:
double dSum = 0;
for(int n = 0; n < 101; n++)
{
dSum += exp((-(xi_n^2*D*bd) / (a^2))* 2* (1 + sin(2*xi_n)/(2*xi_n))^-1 * ((2*PI*q*a)*sin(2*PI*q*a)*cos(xi_n)-xi_n*cos(2*PI*q*a)*sin(xi_n))^2 / ((2*PI*q*a)^2-xi_n^2)^2) + exp((-(zeta_n^2*D*bd) / (a^2))* 2* (1 - sin(2*zeta_n)/(2*zeta_n))^-1 * ((2*PI*q*a)*cos(2*PI*q*a)*sin(zeta_n)-zeta_n*sin(2*PI*q*a)*cos(zeta_n))^2 / ((2*PI*q*a)^2-zeta_n^2)^2) ;
}
E =I0*dSum + K;
I then need to find the transcendental roots in Origin C with the Newton-Raphson algorithm. This looks like:
#include <math.h>
#define MAXIT 100 //max number of iterations
float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2, float xacc)
//using a combination of Newton-Raphson and bisection, find the root of a function bracketed
//between x1 and x2. The root, returned as the function value "rtsafe", will be refined until
//its accuracy is known within +/-xacc. "funcd" is a user-supplied routine that returns both the
//funciton value and the first derivative of the function.
{
void nrerror(char error_text[]);
int j;
float df,dx,dxold,f,fh,fl;
float temp,xh,kl.rts;
(*funcd)(x1,&fl,&df);
(*funcd)(x2,&fh,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
nrerror("Root must be bracketed in rtsafe");
if (fl == 0.0) return x1;
if (fh == 0.0) return x2;
if (fl < 0.0) { //orient the search so that f(x1) < 0
xl=x1;
xh=x2;
} else {
xh=x1;
xl=x2;
}
rts=0.5*(x1+x2); //initialize the guess for root,
dxold=fabs(x2-x1); //the "stepsize before last,"
dx=dxold; //and the last step.
(*funcd)(rts,&f,&df);
for (j=1;j<=MAXIT;j++) { //Loop over allowed iterations.
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
|| (fabs(2.0*f) > fabs(dxold*df))) { //or not decreasing fast enough.
dxold=dx;
dx=0.5*(xh-xl);
rts=xl+dx;
if (xl == rts) return rts; //Change in root is negligible.
} else { //Newton step acceptable. Take it.
dxold=dx;
dx=f/df;
temp=rts;
rts -= dx;
if (temp == rts) return rts;
}
if (fabs(dx) < xacc) return rts; //Convergence criterion.
(*funcd)(rts,&f,&df); //The one new function evaluation per iteration.
if (f < 0.0) //Maintain the bracket on the root.
xl=rts;
else
xh=rts;
}
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0; //Never get here.
}
Continues below! |
| 1 L A T E S T R E P L I E S (Newest First) |
| espenhjo |
Posted - 10/15/2003 : 08:13:21 AM
Is the headerfile <math.h> included in Origin C?
How should I define the functions and the first derivative of the function?
How should it be set up for the incrementing of n to work?
Should I optimize both of the transcendental equations or just one and then calculate the other?
Many big questions! I hope they’re manageable!
Best regards,
Espen
|
|
|
