replaceable function solve
"Solve f_nonlinear(x_zero)=y_zero; f_nonlinear(x_min) - y_zero and f_nonlinear(x_max)-y_zero must have different sign"
import Modelica.Utilities.Streams.error;
extends Modelica.Icons.Function;
input Real y_zero
"Determine x_zero, such that f_nonlinear(x_zero) = y_zero";
input Real x_min
"Minimum value of x";
input Real x_max
"Maximum value of x";
input Real pressure=0.0
"Disregarded variables (here always used for pressure)";
input Real[:] X=
fill(0, 0)
"Disregarded variables (here always used for composition)";
input f_nonlinear_Data f_nonlinear_data
"Additional data for function f_nonlinear";
input Real x_tol=100*Modelica.Constants.eps
"Relative tolerance of the result";
output Real x_zero
"f_nonlinear(x_zero) = y_zero";
protected
constant Real eps=Modelica.Constants.eps
"Machine epsilon";
constant Real x_eps=1e-10
"Slight modification of x_min, x_max, since x_min, x_max are usually exactly at the borders T_min/h_min and then small numeric noise may make the interval invalid";
Real x_min2=x_min - x_eps;
Real x_max2=x_max + x_eps;
Real a=x_min2
"Current best minimum interval value";
Real b=x_max2
"Current best maximum interval value";
Real c
"Intermediate point a <= c <= b";
Real d;
Real e
"b - a";
Real m;
Real s;
Real p;
Real q;
Real r;
Real tol;
Real fa
"= f_nonlinear(a) - y_zero";
Real fb
"= f_nonlinear(b) - y_zero";
Real fc;
Boolean found=false;
algorithm
// Check that f(x_min) and f(x_max) have different sign
fa :=
f_nonlinear(
x_min2,
pressure,
X,
f_nonlinear_data) - y_zero;
fb :=
f_nonlinear(
x_max2,
pressure,
X,
f_nonlinear_data) - y_zero;
fc := fb;
if fa > 0.0
and fb > 0.0
or fa < 0.0
and fb < 0.0
then
error(
"The arguments x_min and x_max to OneNonLinearEquation.solve(..)\n"
+ "do not bracket the root of the single non-linear equation:\n" +
" x_min = " +
String(x_min2) + "\n" + " x_max = " +
String(x_max2)
+ "\n" + " y_zero = " +
String(y_zero) + "\n" +
" fa = f(x_min) - y_zero = " +
String(fa) + "\n" +
" fb = f(x_max) - y_zero = " +
String(fb) + "\n" +
"fa and fb must have opposite sign which is not the case");
end if;
// Initialize variables
c := a;
fc := fa;
e := b - a;
d := e;
// Search loop
while not found
loop
if abs(fc) <
abs(fb)
then
a := b;
b := c;
c := a;
fa := fb;
fb := fc;
fc := fa;
end if;
tol := 2*eps*
abs(b) + x_tol;
m := (c - b)/2;
if abs(m) <= tol
or fb == 0.0
then
// root found (interval is small enough)
found := true;
x_zero := b;
else
// Determine if a bisection is needed
if abs(e) < tol
or abs(fa) <=
abs(fb)
then
e := m;
d := e;
else
s := fb/fa;
if a == c
then
// linear interpolation
p := 2*m*s;
q := 1 - s;
else
// inverse quadratic interpolation
q := fa/fc;
r := fb/fc;
p := s*(2*m*q*(q - r) - (b - a)*(r - 1));
q := (q - 1)*(r - 1)*(s - 1);
end if;
if p > 0
then
q := -q;
else
p := -p;
end if;
s := e;
e := d;
if 2*p < 3*m*q -
abs(tol*q)
and p <
abs(0.5*s*q)
then
// interpolation successful
d := p/q;
else
// use bi-section
e := m;
d := e;
end if;
end if;
// Best guess value is defined as "a"
a := b;
fa := fb;
b := b + (
if abs(d) > tol
then d
else if m > 0
then tol
else -tol);
fb :=
f_nonlinear(
b,
pressure,
X,
f_nonlinear_data) - y_zero;
if fb > 0
and fc > 0
or fb < 0
and fc < 0
then
// initialize variables
c := a;
fc := fa;
e := b - a;
d := e;
end if;
end if;
end while;
end solve;