Hello,
Any body who master a C language to translate math function (non linear equation root) to xojo ?
How many lines of code?
Can you post the code here?
Just post the code so we can have a look.
#include <math.h>
#define MAXIT 100
float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,
float xacc)
{
void nrerror(char error_text);
int j;
float df,dx,dxold,f,fh,fl;
float temp,xh,xl,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) {
xl=x1;
xh=x2;
} else {
xh=x1;
xl=x2;
}
rts=0.5*(x1+x2);
dxold=fabs(x2-x1);
dx=dxold;
(*funcd)(rts,&f,&df);
for (j=1;j<=MAXIT;j++) {
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
|| (fabs(2.0*f) > fabs(dxold*df))) {
dxold=dx;
dx=0.5*(xh-xl);
rts=xl+dx;
if (xl == rts) return rts;
} else {
dxold=dx;
dx=f/df;
temp=rts;
rts -= dx;
if (temp == rts) return rts;
}
if (fabs(dx) < xacc) return rts;
(*funcd)(rts,&f,&df);
if (f < 0.0)
xl=rts;
else
xh=rts;
}
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;
}
#undef MAXIT
sorry , i use a blockquote but something wrong , in the begenning a difference ocurs (code format).
The script is Ok .
remove “;”
remove brackets { and }
put “end if” where they belong
end for where it should be
replace == with =
and it should be good to go almost.
It looks like it’s using the secant method:
I’ve already written Xojo code to do this. It’ll save having to wade through the horror of C code. Too late tonight, but I’ll dig it out tomorrow and post it.
#include <math.h>
#define MAXIT 100
float rtsafe(void (*funcd)(float, float *, float *), float x1, float x2,
float xacc)
{
void nrerror(char error_text[]);
int j;
float df,dx,dxold,f,fh,fl;
float temp,xh,xl,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) {
xl=x1;
xh=x2;
} else {
xh=x1;
xl=x2;
}
rts=0.5*(x1+x2);
dxold=fabs(x2-x1);
dx=dxold;
(*funcd)(rts,&f,&df);
for (j=1;j<=MAXIT;j++) {
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
|| (fabs(2.0*f) > fabs(dxold*df))) {
dxold=dx;
dx=0.5*(xh-xl);
rts=xl+dx;
if (xl == rts) return rts;
} else {
dxold=dx;
dx=f/df;
temp=rts;
rts -= dx;
if (temp == rts) return rts;
}
if (fabs(dx) < xacc) return rts;
(*funcd)(rts,&f,&df);
if (f < 0.0)
xl=rts;
else
xh=rts;
}
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;
}
#undef MAXIT
I said earlier that it appears to use the secant method, but I now think it uses the bisection method, which is similar but less flexible than the secant method. In the bisection method, the two input values must bracket the root; otherwise it won’t converge. Here is my code for the secant method.
Function SecantRoot(xA As double, xB As double, errTolerance As double) as double
'Secant method for finding root of function
'Input parameters are two x values which
'preferably bracket the root value, and if not, are at least near the root
'
'SecantTestFn is the name of the function whose root is to be found
'
const maxit = 100
dim x0 As Double = xA
dim x1 As Double = xB
if SecantTestFn(x0)=0 then return x0
if SecantTestFn(x1)=0 then return x1
dim xn As Double
for i as Integer = 1 to maxit
xn = x1 - SecantTestFn(x1)*(x1-x0)/(SecantTestFn(x1)-SecantTestFn(x0))
if SecantTestFn(xn)=0 or abs(1-x1/xn)<errTolerance then return xn
x0=x1
x1=xn
next
'Failed to find root within specified tolerance. So, return closest value
return xn
End FunctionMuch of the confusing clutter in the posted C code is because one of the calling parameters is a pointer to the function to be minimized. My code has the function name hard coded. This could be changed by passing the function as a delegate.
1 Like
RTSafe is bisection + Newton-Raphson
See pdf page 390 (printed book 366): https://www.cec.uchile.cl/cinetica/pcordero/MC_libros/NumericalRecipesinC.pdf
Yes,
Yes the function posted rtsafe is from NumericalRecipesC.
the algo use the rtsafe to find the root, the author recommand rtsafe instead of NewtonRaphson
.
In the algo the author use a Lagrange Muliplier (optimization) wich is the smalest root of an a 6th deg equation.
I think rtsafe is used to insure the convergeance if a local extremum appears.
gooood , thank’s.
Ah, Numerical Recipes, of course! No wonder the code is unreadable.
rtsafe converted to Xojo. Not tested, and probably won’t work without some modifications, but this should give a rough idea. I have almost no experience with delegates. The last time I used them was many years ago. So hopefully someone will jump in here and correct what I’ve done.
function rtsafe (funcd as delegate, x1 as double, x2 as double) as double
dim j as integer
dim df,dx,dxold,f,fh,fl as double
dim temp,xh,xl, rts as double
funcd.invoke(x1,fl,df)
funcd.invoke(x2,fh,df)
if ((fl > 0.0 and fh > 0.0) or (fl < 0.0 and fh < 0.0)) then
nrerror("Root must be bracketed in rtsafe")
end if
if (fl = 0.0) then return x1
if (fh = 0.0) then return x2
if (fl < 0.0) then
xl = x1
xh = x2
else
xh = x1
xl = x2
end if
rts = 0.5 *(x1 + x2)
dxold = abs(x2 - x1)
dx = dxold
funcd.invoke(rts,f,df)
for j = 1 to MAXIT
if ((((rts - xh) *df - f) *((rts - xl) *df - f) > 0.0) or (fabs(2.0 *f) > fabs(dxold *df))) then
dxold = dx
dx = 0.5 *(xh - xl)
rts = xl + dx
if (xl = rts) then return rts
else
dxold = dx
dx = f / df
temp = rts
rts = rts - dx
if (temp = rts) then return rts
end if
if (fabs(dx) < xacc) then return rts
funcd.invoke(rts, f, df)
if (f < 0.0) then
xl = rts
else
xh = rts
end if
next
end functionHere is a more complete set of the Robert’s conversion.
I felt I needed to make clear the necessity of the ByRefs or it would fail miserably case the user misses it.
Public Sub funcdDelegate(value As Double, Byref result As double, Byref derivative As Double)
Public Function rtsafe(funcDeriv As funcdDelegate, x1 As Double, x2 As Double, maxIter As Integer = 100, accuracy As Double = 3E-8) as Double
Var j As Integer
Var df,dx,dxold,f,fh,fl As Double
Var temp,xh,xl, rts As Double
funcDeriv.Invoke(x1,fl,df)
funcDeriv.Invoke(x2,fh,df)
If ((fl > 0.0 And fh > 0.0) Or (fl < 0.0 And fh < 0.0)) Then
Var ex As New RuntimeException
ex.Message = "Root must be bracketed in rtsafe"
Raise ex
Return 0.0
End If
If (fl = 0.0) Then Return x1
If (fh = 0.0) Then Return x2
If (fl < 0.0) Then
xl = x1
xh = x2
Else
xh = x1
xl = x2
End If
rts = 0.5 *(x1 + x2)
dxold = abs(x2 - x1)
dx = dxold
funcDeriv.Invoke(rts,f,df)
For j = 1 To maxIter
If (((rts - xh) *df - f) *((rts - xl) *df - f) > 0.0) Or (abs(2.0 *f) > abs(dxold *df)) Then
dxold = dx
dx = 0.5 *(xh - xl)
rts = xl + dx
If (xl = rts) Then Return rts
Else
dxold = dx
dx = f / df
temp = rts
rts = rts - dx
If (temp = rts) Then Return rts
End If
If (abs(dx) < accuracy) Then Return rts
funcDeriv.Invoke(rts, f, df)
If (f < 0.0) Then
xl = rts
Else
xh = rts
End If
Next
Var ex As New RuntimeException
ex.Message = "Maximum number of iterations exceeded in rtsafe"
Raise ex
Return 0.0
End Function
Public Sub funcd(value As Double, Byref result As double, Byref derivative As Double)
// funcd example matching the type funcdDelegate
// Do your math here receiving "value" and returning the result and the first derivative
End Sub
Public Sub main()
// Example passing the delegated function funcd above
Var root as double = rtsafe((New funcdDelegate(AddressOf funcd)), x1, x2)
End Sub
2 Likes
Rober Weaver,
the most work is done.
many thank’s.
Thanks. I was going to mention the ByRef requirement, but wasn’t sure of the syntax for the delegate definition. So, I thought it best to leave well enough alone, and let someone else do it properly.
Also, I see you caught the fabs() and changed it to abs(). I just noticed that I’d missed it.
1 Like
Hello Rick,
that what i was trying to do …
thank for you help.
1 Like
Hi Rick,
example to use rtsafe().
// Function F(x)
Fx(X as double) as double
return cos(x)-x * x * x //NewtonRaphson failed to converge to the root with
guess = 0
//First derivative
dFdx(X as double) as double
return -sin(x)-3 * x * x
With bisection [0,1] the convergeance done, root : 0.865479.
Could you show me by an example how to apply this example in rtsafe() ?
Not understand funcdDelegate() …