replaceable record f_nonlinear_Data "Data specific for function f_nonlinear" extends Modelica.Icons.Record; end f_nonlinear_Data;
Type | Name | Default | Description |
---|---|---|---|
Real | y_zero | Determine x_zero, such that f_nonlinear(x_zero) = y_zero | |
Real | x_min | Minimum value of x | |
Real | x_max | Maximum value of x | |
Real | pressure | 0.0 | disregaded variables (here always used for pressure) |
Real | X[:] | fill(0, 0) | disregaded variables (here always used for composition) |
f_nonlinear_Data | f_nonlinear_data | Additional data for function f_nonlinear | |
Real | x_tol | 100*Modelica.Constants.eps | Relative tolerance of the result |
Type | Name | Description |
---|---|---|
Real | x_zero | f_nonlinear(x_zero) = y_zero |
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 "disregaded variables (here always used for pressure)"; input Real[:] X = fill(0,0) "disregaded 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;