| T O P I C R E V I E W |
| riclambo@hotmail.com |
Posted - 07/03/2019 : 07:03:39 AM
Origin Ver. and Service Release (Select Help-->About Origin):
Dear Forum,
I'm trying to fit my data to the following integral.
Here, alpha is the dependent variable, ‘rho’(‘x’in the code below) is the integral independent variable. The parameters I need to fit for are λ （‘l’in the code below）and csda （‘r0’in the code below）. rg is the gyroscopic radius, E0/(q*B)，where ‘B' is the fitting independent variable and E0 and q are constants.
I have followed this example:
https://www.originlab.com/doc/Tutorials/Fitting-Integral-ParaLimit-NAG
The numerator can only be evaluated numerically. The chief difficulty is that it contains a removable singularity. When B=0, the numerator and denominator are equal and alpha = 1. Origin may not realize the convergence and so I've tried to work around it with the following code.
#include <origin.h>
// Add your special include files here.
// For example, if you want to fit with functions from the NAG library,
// add the header file for the NAG functions here.
#include <OC_nag.h>
// Add code here for other Origin C functions that you want to define in this file,
// and access in your fitting function.
struct user
{
double l, E0, q, mu, fitX; // fitX the independent variable of fitting function
};
static double NAG_CALL f_callback(double x, Nag_User *comm) // x is the independent variable of the integrand
{
struct user *sp = (struct user *)(comm->p);
double ll, E0E0, qq, mumu, fitX; // temp variable to accept the parameters in the Nag_User communication struct
ll = sp->l;
E0E0 = sp->E0;
qq = sp->q;
mumu = sp->mu;
fitX = sp->fitX;
return x*exp(mumu*x - (2*E0E0/(qq*fitX*ll))*asin(x/(2*E0E0*E0E0/(qq*fitX)))) ;
}
// You can access C functions defined in other files, if those files are loaded and compiled
// in your workspace, and the functions have been prototyped in a header file that you have
// included above.
// You can access NLSF object methods and properties directly in your function code.
// You should follow C-language syntax in defining your function.
// For instance, if your parameter name is P1, you cannot use p1 in your function code.
// When using fractions, remember that integer division such as 1/2 is equal to 0, and not 0.5
// Use 0.5 or 1/2.0 to get the correct value.
// For more information and examples, please refer to the "User-Defined Fitting Function"
// section of the Origin Help file.
//----------------------------------------------------------
//
void _nlsfalphavslowB(
// Fit Parameter(s):
double r0, double l,
// Independent Variable(s):
double B,
// Dependent Variable(s):
double& y)
{
const double E0=10;
const double q=3;
const double mu=0.01579;
// Beginning of editable part
// Beginning of editable part
double epsabs = 0.00001, epsrel = 0.0000001, result, abserr;
Integer max_num_subint = 500;
// you may use epsabs and epsrel and this quantity to enhance your desired precision
// when not enough precision encountered
Nag_QuadProgress qp;
static NagError fail;
// the parameters parameterize the integrand can be input to the call_back function
// through the Nag_User communication struct
Nag_User comm;
struct user s;
s.l = l;
s.E0 = E0;
s.q = q;
s.mu = mu;
s.fitX = B;
comm.p = (Pointer)&s;
if (B==0)
{
(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
}
else{
d01sjc(f_callback, 0, r0, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
}
// you may want to exam the error by printing out error message, just uncomment the following lines
// if (fail.code != NE_NOERROR)
// printf("%s\n", fail.message);
// For the error other than the following three errors which are due to bad input parameters
// or allocation failure NE_INT_ARG_LT NE_BAD_PARAM NE_ALLOC_FAIL
// You will need to free the memory allocation before calling the integration routine again to
// avoid memory leakage
if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
{
NAG_FREE(qp.sub_int_beg_pts);
NAG_FREE(qp.sub_int_end_pts);
NAG_FREE(qp.sub_int_result);
NAG_FREE(qp.sub_int_error);
}
y = result/(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
// End of editable part
}
I have tried to get around it by the use of an if/else condition. But when I test the function I can see it does not give me the right value at B=0. Are there any more experienced users who can help.
The data:
data = {{0, 1}, {0.3, 0.999552770379519}, {0.4,
0.998805851062253}, {0.5, 0.998066565041057}, {0.6,
0.997060919995234}};
Regards,
RL
RL |
| 13 L A T E S T R E P L I E S (Newest First) |
| yuki_wu |
Posted - 07/28/2019 : 10:41:24 PM
Hi Jeff,
This quick help might help:
https://www.originlab.com/doc/Quick-Help/Summation-DoubleIntegral-Func
BTW, could you please open a new topic for different question next time?
Regards,
Yuki
OriginLab
|
| Kavitamarne |
Posted - 07/26/2019 : 04:42:42 AM
Excellent thread! Thank you! |
| riclambo@hotmail.com |
Posted - 07/09/2019 : 05:38:03 AM
I finally got it to work. Here is the code:
#include <origin.h>
#include <origin.h>
// Add your special include files here.
// For example, if you want to fit with functions from the NAG library,
// add the header file for the NAG functions here.
#include <OC_nag.h>
// Add code here for other Origin C functions that you want to define in this file,
// and access in your fitting function.
struct user
{
double l, E0, q, mu, fitX; // fitX the independent variable of fitting function
};
static double NAG_CALL f_callback(double x, Nag_User *comm) // x is the independent variable of the integrand
{
struct user *sp = (struct user *)(comm->p);
double ll, E0E0, qq, mumu, fitX; // temp variable to accept the parameters in the Nag_User communication struct
ll = sp->l;
E0E0 = sp->E0;
qq = sp->q;
mumu = sp->mu;
fitX = sp->fitX;
return x*exp(mumu*x - (2*E0E0/(qq*fitX*ll))*asin(x/(2*E0E0/(qq*fitX)))) ;
}
// You can access C functions defined in other files, if those files are loaded and compiled
// in your workspace, and the functions have been prototyped in a header file that you have
// included above.
// You can access NLSF object methods and properties directly in your function code.
// You should follow C-language syntax in defining your function.
// For instance, if your parameter name is P1, you cannot use p1 in your function code.
// When using fractions, remember that integer division such as 1/2 is equal to 0, and not 0.5
// Use 0.5 or 1/2.0 to get the correct value.
// For more information and examples, please refer to the "User-Defined Fitting Function"
// section of the Origin Help file.
//----------------------------------------------------------
//
void _nlsfalphavslowB_test(
// Fit Parameter(s):
double r0, double l,
// Independent Variable(s):
double B,
// Dependent Variable(s):
double& y)
{
const double E0=10;
const double q=3;
const double mu=0.01579;
// Beginning of editable part
// Beginning of editable part
double epsabs = 0.00001, epsrel = 0.0000001, result, abserr, result1;
Integer max_num_subint = 1000;
// you may use epsabs and epsrel and this quantity to enhance your desired precision
// when not enough precision encountered
Nag_QuadProgress qp;
static NagError fail;
// the parameters parameterize the integrand can be input to the call_back function
// through the Nag_User communication struct
Nag_User comm;
struct user s;
s.l = l;
s.E0 = E0;
s.q = q;
s.mu = mu;
s.fitX = B;
comm.p = (Pointer)&s;
d01sjc(f_callback, 0, r0, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
d01sjc(f_callback1, 0, r0, epsabs, epsrel, max_num_subint, &result1, &abserr, &qp, &comm, &fail);
// you may want to exam the error by printing out error message, just uncomment the following lines
// if (fail.code != NE_NOERROR)
// printf("%s\n", fail.message);
// For the error other than the following three errors which are due to bad input parameters
// or allocation failure NE_INT_ARG_LT NE_BAD_PARAM NE_ALLOC_FAIL
// You will need to free the memory allocation before calling the integration routine again to
// avoid memory leakage
if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
{
NAG_FREE(qp.sub_int_beg_pts);
NAG_FREE(qp.sub_int_end_pts);
NAG_FREE(qp.sub_int_result);
NAG_FREE(qp.sub_int_error);
}
if(B==0) {
y = 1;
} else {
y = result/(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0)))/((-1
+l*mu)^2));
}
// End of editable part
}
Thanks Castiel. Thanks Yuki. It was just a question of putting the if/else clause in the right place.
RL |
| riclambo@hotmail.com |
Posted - 07/09/2019 : 12:03:27 AM
I more or less understand the problem with the code now. But what is the solution？ How can I ‘force’y=1 at B=0? What about a piece fit with two curves？
RL |
| Castiel |
Posted - 07/08/2019 : 11:48:46 PM
quote:
Originally posted by riclambo@hotmail.com
I really appreciate your taking a look, Castiel, since this is a struggle for a weak programmer like myself. As you say, the following lines seem to be the problem:
if (B==0)
{
(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
}
else{
d01sjc(f_callback, 0, r0, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
}
However, what I'm trying to do is make "y=1" when B=0. If you look down to near the end of the code, I've written:
y = result/(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
// End of editable part
In other words, I'd hoped numerator and denominator would cancel to give y=1 at B=0. This was the idea, however, I'm not quite able to execute it.
RL
Given B = 0, the following code does nothing:
if (B==0)
{
(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
}
"result" is left unchanged as 0, therefore y is either 0 or NaN (depends on other parameters values) but not 1.
&&&&&&&&&
&&&
&&
& _____ ___________
II__|[] | | I I |
| |_|_ I I _|
< OO----OOO OO---OO
**********************************************************
|
| riclambo@hotmail.com |
Posted - 07/08/2019 : 10:35:39 PM
I really appreciate your taking a look, Castiel, since this is a struggle for a weak programmer like myself. As you say, the following lines seem to be the problem:
if (B==0)
{
(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
}
else{
d01sjc(f_callback, 0, r0, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
}
However, what I'm trying to do is make "y=1" when B=0. If you look down to near the end of the code, I've written:
y = result/(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
// End of editable part
In other words, I'd hoped numerator and denominator would cancel to give y=1 at B=0. This was the idea, however, I'm not quite able to execute it.
RL |
| Castiel |
Posted - 07/08/2019 : 10:14:24 PM
quote:
Originally posted by riclambo@hotmail.com
Got it,Castiel. However, I was just using sin(x)/x to illustrate my main point, which is how to fit with a function that has a removable singularity. The function that I am genuinely interested in is the one that appears in my first post.
RL
"When B=0, the numerator and denominator are equal and alpha = 1. Origin may not realize the convergence and so I've tried to work around it with the following code."
if (B==0)
{
(l*(l + ((-r0 +l*(-1 + mu*r0))*exp((mu - 1/l)*r0))/((-1 +l*mu)^2)));
}
That seems to be a problem.
If the result is incorrect, debug your code.
&&&&&&&&&
&&&
&&
& _____ ___________
II__|[] | | I I |
| |_|_ I I _|
< OO----OOO OO---OO
**********************************************************
|
| riclambo@hotmail.com |
Posted - 07/08/2019 : 9:38:39 PM
Got it,Castiel. However, I was just using sin(x)/x to illustrate my main point, which is how to fit with a function that has a removable singularity. The function that I am genuinely interested in is the one that appears in my first post.
RL |
| Castiel |
Posted - 07/08/2019 : 11:07:03 AM
quote:
Originally posted by riclambo@hotmail.com
Hi Yuki,
Your reply did provide some help. As you said, when x= 0, sin[x]/x will be a missing value. If I try to integrate it according to your code, the integral will also have a missing value.
However, I can input the value x=0, sin[x]/x if I use this code:
newbook;
col(a)={-5:0.5:5};
col(b)=sin(col(a))/col(a);
cell(11, 2) = 1;
col(c)=IntegrateXY(col(a), col(b));
col(c)[C]$ = "Area";
Please try it to see the difference. It will give me a continous function and a continuous integral. This is what I'm trying to do with my integral. I'm trying to find out how to make the program realise that when rho=0, alpha =1.
I submitted this question to tech@originlab.com, but I have not recieved a reply yet.
Probably what I need to do is a piecewise fit. Something like:
https://www.originlab.com/doc/Tutorials/Fitting-Piecewise-Linear
RL
That's a special function.
https://www.originlab.com/doc/LabTalk/ref/Sin-integral-func
&&&&&&&&&
&&&
&&
& _____ ___________
II__|[] | | I I |
| |_|_ I I _|
< OO----OOO OO---OO
**********************************************************
|
| riclambo@hotmail.com |
Posted - 07/08/2019 : 07:13:56 AM
Hi Yuki,
Your reply did provide some help. As you said, when x= 0, sin[x]/x will be a missing value. If I try to integrate it according to your code, the integral will also have a missing value.
However, I can input the value x=0, sin[x]/x if I use this code:
newbook;
col(a)={-5:0.5:5};
col(b)=sin(col(a))/col(a);
cell(11, 2) = 1;
col(c)=IntegrateXY(col(a), col(b));
col(c)[C]$ = "Area";
Please try it to see the difference. It will give me a continous function and a continuous integral. This is what I'm trying to do with my integral. I'm trying to find out how to make the program realise that when rho=0, alpha =1.
I submitted this question to tech@originlab.com, but I have not recieved a reply yet.
Probably what I need to do is a piecewise fit. Something like:
https://www.originlab.com/doc/Tutorials/Fitting-Piecewise-Linear
RL |
| yuki_wu |
Posted - 07/07/2019 : 11:20:15 PM
Hi,
When x = 0, sin(x)/x shoule be a missing value. Missing value will be ignored when we calculate the integration.
For example, you could execute the following LabTalk script to check the result:
newbook;
col(a)={-5:0.5:5};
col(b)=sin(col(a))/col(a);
col(c)=IntegrateXY(col(a), col(b));
col(c)[C]$ = "Area";
Hope it can be some help.
Regards,
Yuki
OriginLab
|
| riclambo@hotmail.com |
Posted - 07/04/2019 : 10:29:20 AM
Thanks for replying so quickly. I will send you the materials you requested, but I would like to give my basic question one more try.
Boiling it down, it comes to something like this. How would you fit the following function numerically to your data?
It looks like there is a singularity at x=0, but in reality it is a removable singularity. At that point the integrand is equal to 1.
This is the same problem I'm facing. At rho=0, the value of alpha should be 1.
RL |
| yuki_wu |
Posted - 07/04/2019 : 01:42:29 AM
Hi,
Sorry, it is not easy to understand your question and figure out what went wrong from the long description above.
You better send us your FDF file and data via tech@originlab.com. Also tell us your serial number and which version you are using.
Regards,
Yuki
OriginLab
|
|
|
