Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module num_lib
21 : ! various numerical routines
22 :
23 : use const_def, only: dp, qp, iln10
24 : use math_lib
25 : use num_def
26 : use mod_integrate
27 :
28 : ! NOTE: because of copyright restrictions,
29 : ! mesa doesn't use any routines from Numerical Recipes.
30 :
31 : implicit none
32 :
33 : private
34 : public :: find0
35 : public :: find0_quadratic
36 : public :: binary_search
37 : public :: binary_search_sg
38 : public :: safe_root
39 : public :: safe_root_with_guess
40 : public :: safe_root_with_brackets
41 : public :: safe_root_without_brackets
42 : public :: null_mas
43 : public :: null_fcn
44 : public :: null_jac
45 : public :: null_sjac
46 : public :: null_fcn_blk_dble
47 : public :: null_jac_blk_dble
48 : public :: brent_safe_zero
49 : public :: brent_local_min
50 : public :: brent_global_min
51 : public :: default_bdomain
52 : public :: default_force_another_iter
53 : public :: default_inspectb
54 : public :: default_set_primaries
55 : public :: default_set_secondaries
56 : public :: default_size_equ
57 : public :: default_sizeb
58 : public :: default_xdomain
59 : public :: integrate
60 : public :: bobyqa
61 : public :: dop853
62 : public :: dop853_work_sizes
63 : public :: isolve
64 : public :: isolve_work_sizes
65 : public :: dopri5
66 : public :: dopri5_work_sizes
67 : public :: newton
68 : public :: newton_work_sizes
69 : public :: newuoa
70 : public :: find_max_quadratic
71 : public :: look_for_brackets
72 : public :: nm_simplex
73 : public :: qsort
74 : public :: qsort_string_index
75 : public :: iln10
76 : public :: dfridr
77 : public :: solver_option
78 : public :: linear_interp
79 : public :: two_piece_linear_coeffs
80 : public :: simplex_info_str
81 : public :: simplex_op_code
82 :
83 : contains
84 :
85 : ! tests of numerical derivative
86 : include "num_dfridr.dek"
87 :
88 : ! safe root finding
89 : ! uses alternating bisection and inverse parabolic interpolation
90 : ! also have option to use derivative as accelerator (newton method)
91 : include "num_safe_root.dek"
92 :
93 : ! solvers for ODEs and DAEs.
94 :
95 : ! sometimes you just want a simple runge-kutta
96 : include "num_rk2.dek"
97 : include "num_rk4.dek"
98 :
99 : ! but there are lots of fancier options too.
100 :
101 : ! selections from the Hairer family of ODE/DAE integrators.
102 : ! from Ernst Hairer's website: http://www.unige.ch/~hairer/
103 :
104 : ! explicit ODE solvers based on methods of Dormand and Prince (for non-stiff problems)
105 :
106 : ! explicit Runge-Kutta ODE integrator of order 5
107 : ! with dense output of order 4
108 : include "num_dopri5.dek" ! "DOramand PRInce order 5"
109 :
110 : ! explicit Runge-Kutta ODE integrator of order 8
111 : ! with dense output of order 7
112 : include "num_dop853.dek" ! "DOramand Prince order 8(5, 3)"
113 :
114 : ! both integrators have automatic step size control and monitoring for stiffness.
115 :
116 : ! For a description see:
117 : ! Hairer, Norsett and Wanner (1993):
118 : ! Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition.
119 : ! Springer Series in Comput. Math., vol. 8.
120 : ! http://www.unige.ch/~hairer/books.html
121 :
122 : ! implicit solvers (for stiff problems)
123 :
124 : ! there are a bunch of implicit solvers to pick from (listed below),
125 : ! but they all have pretty much the same arguments,
126 : ! so I've provided a general routine, called "isolve", that let's you
127 : ! pass an extra argument to specify which one of the particular solvers
128 : ! that you want to use.
129 : include "num_isolve.dek"
130 :
131 : ! if possible, you should write your code to call isolve
132 : ! rather than calling one of the particular solvers.
133 : ! only call a specific solver if you need a feature it provides
134 : ! that isn't supported by isolve.
135 :
136 : ! you can find an example program using isolve in num/test/src/sample_ode_solver.f
137 :
138 : ! the implicit solver routines
139 :
140 : ! for detailed descriptions of these routines see:
141 : ! Hairer and Wanner (1996):
142 : ! Solving Ordinary Differential Equations.
143 : ! Stiff and Differential-Algebraic Problems. 2nd edition.
144 : ! Springer Series in Comput. Math., vol. 14.
145 : ! http://www.unige.ch/~hairer/books.html
146 :
147 : ! linearly implicit Runge-Kutta method (Rosenbrock)
148 : include "num_ros2.dek" ! L-stable; 2 stages; order 2, 2 function evaluations.
149 : ! ros2 is suitable for use with approximate jacobians such as from numerical differences.
150 : ! ros2 is designed for use with Strang splitting and is reported to be able to cope
151 : ! with large time steps and artificial transients introduced at the beginning of split intervals.
152 : ! see Verwer et al, "A Second-Order Rosenbrock Method Applied to Photochemical Dispersion Problems",
153 : ! SIAM J. Sci. Comput. (20), 1999, 1456-1480.
154 :
155 : include "num_rose2.dek" ! L-stable; 3 stages; order 2, 3 function evaluations.
156 : ! rose2 is unique among the implicit solvers in that the final function evaluation
157 : ! uses the solution vector for the step.
158 :
159 : include "num_rodas3.dek" ! L-stable; 4 stages; order 3, 3 function evaluations.
160 : include "num_rodas4.dek" ! L-stable; 6 stages; order 4, 6 function evaluations.
161 :
162 : ! 3rd order; for parabolic equations.
163 : include "num_ros3p.dek" ! A-stable; 3 stages; order 3, 2 function evaluations.
164 : include "num_ros3pl.dek" ! L-stable; 4 stages; order 3, 3 function evaluations.
165 : ! 4th order; for parabolic equations.
166 : include "num_rodasp.dek" ! L-stable; 6 stages; order 4, 6 function evaluations.
167 :
168 : include "num_solvers_options.dek"
169 :
170 : ! which implicit solver should you use?
171 :
172 : ! somewhat surprisingly, in some cases the solvers
173 : ! that work well at high tolerances will fail with low
174 : ! tolerances and vice-versa. so you need to match
175 : ! the solver to the problem.
176 :
177 : ! your best bet is to try them all on some typical cases.
178 : ! happily this isn't too hard to do since they all
179 : ! use the same function arguments and have (almost)
180 : ! identical calling sequences and options.
181 :
182 : ! flexible choice of linear algebra routines
183 :
184 : ! the solvers need to solve linear systems.
185 : ! this is typically done by first factoring the matrix A
186 : ! and then repeatedly using the factored form to solve
187 : ! A*x=b for various vectors b.
188 :
189 : ! rather than build in a particular matrix solver,
190 : ! the mesa versions of the solvers take as arguments
191 : ! routines to perform these tasks. the mesa/mtx package
192 : ! includes several choices for implementations of the
193 : ! required routines.
194 :
195 : ! dense, banded, or sparse matrix
196 :
197 : ! All the packages allow the matrix to be in dense or banded form.
198 : ! the choice of sparse matrix package is not fixed by the solvers.
199 : ! the only constraint is the the sparse format must be either
200 : ! compressed row sparse or compressed column sparse.
201 : ! the mesa/mtx package comes with one option for a sparse
202 : ! package (based on SPARSKIT), and also has hooks for another (Super_LU).
203 : ! Since the sparse routines are passed as arguments to the solvers,
204 : ! it is possible to experiment with different linear algebra
205 : ! packages without a great deal of effort.
206 :
207 : ! analytical or numerical jacobian
208 :
209 : ! to solve M*y' = f(y), the solvers need to have the jacobian matrix, df/dy.
210 : ! the jacobian can either be calculated analytically by a user supplied routine,
211 : ! or the solver can form a numerical difference estimate by repeatedly
212 : ! evaluating f(y) with slightly different y's. Such numerical jacobians
213 : ! are supported by all the solvers for both dense and banded matrix forms.
214 : ! For the sparse matrix case, only analytical jacobians are allowed.
215 :
216 : ! NOTE: for most implicit solvers, the accuracy of the jacobian influences
217 : ! the rate of convergence, but doesn't impact the accuracy of the solution.
218 : ! however, the rodas solvers are an exception to this rule.
219 : ! they are based on the rosenbrock method which replaces the newton iteration
220 : ! by formulas that directly use the jacobian in the formula for
221 : ! the result. as a result, the rodas solvers depend on having
222 : ! accurate jacobians in order to produce accurate results.
223 :
224 : ! explicit or implicit ODE systems
225 :
226 : ! systems of the form y' = f(y) are called "explicit ODE systems".
227 :
228 : ! systems of the form M*y' = f(y), with M not equal to the identity matrix,
229 : ! are called "implicit ODE systems".
230 :
231 : ! in addition to the usual explicit systems,
232 : ! the solvers can also handle implicit ODE systems
233 : ! in which M is an arbitrary constant matrix,
234 : ! even including the case of M singular.
235 :
236 : ! for M non-constant, see the discussion of "problems with special structure"
237 :
238 : ! problems with special structure
239 :
240 : ! 3 special cases can be handled easily
241 :
242 : ! case 1, second derivatives: y'' = f(t, y, y')
243 : ! case 2, nonconstant matrix: C(x, y)*y' = f(t, y)
244 : ! case 3, both of the above: C(x, y)*y'' = f(t, y, y')
245 :
246 : ! these all work by adding auxiliary variables to the problem and
247 : ! converting back to the standard form with a constant matrix M.
248 :
249 : ! case 1: y'' = f(t, y, y')
250 : ! after add auxiliary variables z, this becomes
251 : ! y' = z
252 : ! z' = f(t, y, z)
253 :
254 : ! case 2: C(x, y)*y' = f(t, y)
255 : ! after add auxiliary variables z, this becomes
256 : ! y' = z
257 : ! 0 = C(x, y)*z - f(t, y)
258 :
259 : ! case 3: C(x, y)*y'' = f(t, y, y')
260 : ! after add auxiliary variables z and u, this becomes
261 : ! y' = z
262 : ! z' = u
263 : ! 0 = C(x, y)*u - f(t, y, z)
264 :
265 : ! The last two cases take advantage of the ability to have M singular.
266 :
267 : ! If the matrix for df/dy is dense in these special cases, all the solvers
268 : ! can reduce the cost of the linear algebra operations by special treatment
269 : ! of the auxiliary variables.
270 :
271 : ! "projection" of solution to valid range of values.
272 :
273 : ! it is often the case that the n-dimensional solution
274 : ! is actually constrained to a subspace of full n dimensional
275 : ! space of numbers. The proposed solutions at each step
276 : ! need to be projected back to the allowed subspace in order
277 : ! to maintain valid results. The routines all provide for this
278 : ! option by calling a "solout" routine, supplied by the user,
279 : ! after every accepted step. The user's solout routine can modify
280 : ! the solution y before returning to continue the integration.
281 :
282 : ! "dense output"
283 :
284 : ! the routines provide estimates of the solution over entire step.
285 : ! useful for tabulating the solution at prescribed points
286 : ! or for smooth graphical presentation of the solution.
287 : ! also very useful for "event location" -- e.g., at what
288 : ! value of x do we get a solution y(x) s.t. some relation
289 : ! g(x, y(x))=0. The dense output option is very helpful here.
290 : ! All of the solvers support dense output.
291 : ! BTW: there is typically a certain overhead associated with
292 : ! providing the option for dense output, so don't request it
293 : ! unless you'll really be using it.
294 :
295 : ! here is a special version of safe_root for use with dense output "solout" routines
296 : include "num_solout_root.dek"
297 :
298 : ! "null" implementations of routines used by the solvers
299 : ! are for use in cases in which you aren't actually using the routine,
300 : ! but something must be provided for the required argument.
301 :
302 0 : subroutine null_fcn(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
303 : integer, intent(in) :: n, lrpar, lipar
304 : real(dp), intent(in) :: x, h
305 : real(dp), intent(inout) :: y(:) ! (n)
306 : ! okay to edit y if necessary (e.g., replace negative values by zeros)
307 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
308 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
309 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
310 : integer, intent(out) :: ierr ! nonzero means retry with smaller timestep.
311 0 : f = 0; ierr = 0
312 0 : end subroutine null_fcn
313 :
314 0 : subroutine null_fcn_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
315 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
316 : real(dp), intent(in) :: x, h
317 : real(dp), intent(inout), pointer :: y(:)
318 : ! (n) okay to edit y if necessary (e.g., replace negative values by zeros)
319 : real(dp), intent(inout), pointer :: f(:) ! (n) dy/dx
320 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
321 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
322 : integer, intent(out) :: ierr ! nonzero means retry with smaller timestep.
323 0 : f = 0; ierr = 0
324 0 : end subroutine null_fcn_blk_dble
325 :
326 0 : subroutine null_jac(n, x, h, y, f, dfy, ldfy, lrpar, rpar, lipar, ipar, ierr)
327 : integer, intent(in) :: n, ldfy, lrpar, lipar
328 : real(dp), intent(in) :: x, h
329 : real(dp), intent(inout) :: y(:) ! (n)
330 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
331 : real(dp), intent(inout) :: dfy(:, :) ! (ldfy, n)
332 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
333 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
334 : integer, intent(out) :: ierr ! nonzero means terminate integration
335 0 : f = 0; dfy = 0; ierr = 0
336 0 : end subroutine null_jac
337 :
338 0 : subroutine null_jac_blk_dble(n, caller_id, nvar, nz, x, h, y, f, lblk, dblk, ublk, lrpar, rpar, lipar, ipar, ierr)
339 : integer, intent(in) :: n, caller_id, nvar, nz, lrpar, lipar
340 : real(dp), intent(in) :: x, h
341 : real(dp), intent(inout), pointer :: y(:) ! (n)
342 : real(dp), intent(inout), pointer :: f(:) ! (n) dy/dx
343 : real(dp), dimension(:), pointer, intent(inout) :: lblk, dblk, ublk ! =(nvar,nvar,nz)
344 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
345 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
346 : integer, intent(out) :: ierr ! nonzero means terminate integration
347 0 : f = 0; y = 0; ierr = 0
348 0 : end subroutine null_jac_blk_dble
349 :
350 0 : subroutine null_sjac(n, x, h, y, f, nzmax, ia, ja, values, lrpar, rpar, lipar, ipar, ierr)
351 : ! sparse jacobian. format either compressed row or compressed column.
352 : integer, intent(in) :: n, nzmax, lrpar, lipar
353 : real(dp), intent(in) :: x, h
354 : real(dp), intent(inout) :: y(:) ! (n)
355 : real(dp), intent(inout) :: f(:) ! (n) ! dy/dx
356 : integer, intent(inout) :: ia(:) ! (n+1)
357 : integer, intent(inout) :: ja(:) ! (nzmax)
358 : real(dp), intent(inout) :: values(:) ! (nzmax)
359 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
360 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
361 : integer, intent(out) :: ierr ! nonzero means terminate integration
362 0 : f = 0; values = 0; ia = 0; ja = 0; ierr = 0
363 0 : end subroutine null_sjac
364 :
365 0 : subroutine null_mas(n, am, lmas, lrpar, rpar, lipar, ipar)
366 : integer, intent(in) :: n, lmas, lrpar, lipar
367 : real(dp), intent(inout) :: am(:, :) ! (lmas, n)
368 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
369 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
370 0 : am = 0
371 0 : end subroutine null_mas
372 :
373 : subroutine null_solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
374 : integer, intent(in) :: nr, n, lrpar, lipar
375 : real(dp), intent(in) :: xold, x
376 : real(dp), intent(inout) :: y(:) ! (n)
377 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
378 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
379 : real(dp), intent(inout), target :: rwork(*)
380 : integer, intent(inout), target :: iwork(*)
381 : interface
382 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
383 : use const_def, only: dp
384 : implicit none
385 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
386 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
387 : real(dp), intent(inout), target :: rwork(*)
388 : integer, intent(inout), target :: iwork(*)
389 : integer, intent(out) :: ierr
390 : end function interp_y
391 : end interface
392 : integer, intent(out) :: irtrn
393 : irtrn = 0
394 : end subroutine null_solout
395 :
396 0 : subroutine null_dfx(n, x, y, fx, lrpar, rpar, lipar, ipar, ierr)
397 : integer, intent(in) :: n, lrpar, lipar
398 : real(dp), intent(in) :: x, y(:) ! (n)
399 : real(dp), intent(inout) :: fx(:) ! (n)
400 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
401 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
402 : integer, intent(out) :: ierr
403 0 : ierr = 0
404 0 : fx = 0
405 0 : end subroutine null_dfx
406 :
407 : ! Newton-Raphson iterative solver for nonlinear systems
408 : ! square or banded
409 : ! analytic or numerical difference jacobian
410 : ! where possible, reuses jacobian to improve efficiency
411 : ! uses line search method to improve "global" convergence
412 : include "num_newton.dek"
413 :
414 : ! minimize scalar function of many variables without using derivatives.
415 :
416 : ! NEWUOA
417 : ! This subroutine seeks the least value of a function of many variables,
418 : ! by applying a trust region method that forms quadratic models by
419 : ! interpolation. There is usually some freedom in the interpolation
420 : ! conditions, which is taken up by minimizing the Frobenius norm of
421 : ! the change to the second derivative of the model, beginning with the
422 : ! zero matrix.
423 : ! by M.J.D. Powell (mjdp@cam.ac.uk)
424 : ! M.J.D. Powell, "Developments of NEWUOA for unconstrained minimization without derivatives",
425 : ! Department of Applied Mathematics and Theoretical Physics, Cambridge, England, report NA05, 2007.
426 : include "num_newuoa.dek"
427 :
428 : ! BOBYQA
429 : ! Similar to NEWUOA, but the values of the variables are constrained
430 : ! by upper and lower bounds.
431 : ! by M.J.D. Powell (mjdp@cam.ac.uk)
432 : ! M.J.D. Powell, "The BOBYQA algorithm for bound constrained optimization without derivatives",
433 : ! Department of Applied Mathematics and Theoretical Physics, Cambridge, England, report NA06, 2009.
434 : include "num_bobyqa.dek"
435 :
436 : ! Nelder-Mead Simplex Method
437 : ! doesn't use interpolation, so robust with noise data.
438 : include "num_simplex.dek"
439 : ! Nelder, J. A. and Mead, R.
440 : ! "A Simplex Method for Function Minimization."
441 : ! Comput. J. 7, 308-313, 1965.
442 :
443 : ! global or local minimum of scalar function of 1 variable
444 : include "num_brent.dek"
445 :
446 : ! QuickSort. ACM Algorithm 402, van Emden, 1970
447 : ! mesa's implementation from Joseph M. Krahn
448 : ! http://fortranwiki.org/fortran/show/qsort_inline
449 :
450 3 : subroutine qsort(index, n, vals)
451 : use mod_qsort, only: sortp_dp
452 : integer :: index(:), n
453 : real(dp) :: vals(:)
454 3 : call sortp_dp(n, index, vals)
455 3 : end subroutine qsort
456 :
457 : subroutine qsort_strings(index, n, strings)
458 : use mod_qsort, only: sortp_string
459 : integer :: index(:), n
460 : character(len=*), intent(in) :: strings(:)
461 : call sortp_string(n, index, strings)
462 : end subroutine qsort_strings
463 :
464 19 : subroutine qsort_string_index(index, n, string_index, strings)
465 : use mod_qsort, only: sortp_string_index
466 : integer :: index(:), n
467 : integer, intent(in) :: string_index(:) ! (n)
468 : character(len=*), intent(in) :: strings(:) ! 1..maxval(string_index)
469 19 : call sortp_string_index(n, index, string_index, strings)
470 19 : end subroutine qsort_string_index
471 :
472 : ! random numbers
473 : real(dp) function get_dp_uniform_01(seed)
474 : ! returns a unit pseudorandom real(dp)
475 : use mod_random, only: r8_uniform_01
476 : integer(kind=4) :: seed
477 : get_dp_uniform_01 = r8_uniform_01(seed)
478 : end function get_dp_uniform_01
479 :
480 : function get_i4_uniform(a, b, seed)
481 : ! The pseudorandom integer will be scaled to be uniformly distributed
482 : ! between a and b.
483 : use mod_random, only: i4_uniform
484 : integer(kind=4) :: a, b, seed, get_i4_uniform
485 : get_i4_uniform = i4_uniform(a, b, seed)
486 : end function get_i4_uniform
487 :
488 : subroutine get_perm_uniform(n, base, seed, p)
489 : ! selects a random permutation of n integers
490 : use mod_random, only: perm_uniform
491 : integer(kind=4) :: n
492 : integer(kind=4) :: base
493 : integer(kind=4) :: p(n)
494 : integer(kind=4) :: seed
495 : call perm_uniform(n, base, seed, p)
496 : end subroutine get_perm_uniform
497 :
498 : subroutine get_seed_for_random(seed)
499 : ! returns a seed for the random number generator
500 : use mod_random, only: get_seed
501 : integer(kind=4) :: seed
502 : call get_seed(seed)
503 : end subroutine get_seed_for_random
504 :
505 : ! binary search
506 : include "num_binary_search.dek"
507 :
508 0 : real(dp) function linear_interp(x1, y1, x2, y2, x)
509 : real(dp), intent(in) :: x1, y1, x2, y2, x
510 0 : if (x2 == x1) then
511 0 : linear_interp = (y1 + y2)/2
512 : else
513 0 : linear_interp = y1 + (y2 - y1)*(x - x1)/(x2 - x1)
514 : end if
515 0 : end function linear_interp
516 :
517 16880 : real(dp) function find0(xx1, yy1, xx2, yy2) result(x)
518 : ! find x between xx1 and xx2 s.t. linear_interp(xx1, yy1, xx2, yy2, x) == 0
519 : real(dp), intent(in) :: xx1, yy1, xx2, yy2
520 16880 : real(dp) :: a, b
521 16880 : a = (xx1*yy2) - (xx2*yy1)
522 16880 : b = yy2 - yy1
523 16880 : if (((abs(a) >= abs(b)*1d99) .and. ((yy1 >= 0d0 .and. yy2 <= 0d0) .or. (yy1 <= 0d0 .and. yy2 >= 0d0)))) then
524 0 : x = 0.5d0*(xx1 + xx2)
525 16880 : else if (b == 0d0) then
526 0 : x = 0.5d0*(xx1 + xx2)
527 : else
528 16880 : x = a/b
529 : end if
530 16880 : if (yy1*yy2 <= 0) then ! sanity check
531 16650 : if (x > max(xx1, xx2)) x = max(xx1, xx2)
532 16650 : if (x < min(xx1, xx2)) x = min(xx1, xx2)
533 : end if
534 16880 : end function find0
535 :
536 2 : real(dp) function find0_quadratic(xx1, yy1, xx2, yy2, xx3, yy3, ierr) result(x)
537 : ! find x between xx1 and xx3 s.t. quad_interp(xx1, yy1, xx2, yy2, xx3, yy3, x) == 0
538 : ! xx2 between xx1 and xx3; yy1 and yy3 different sign; yy2 between yy1 and yy3.
539 : real(dp), intent(in) :: xx1, yy1, xx2, yy2, xx3, yy3
540 : integer, intent(out) :: ierr
541 2 : real(dp) :: a, b, s2, denom
542 2 : ierr = 0; x = 0
543 : s2 = (xx3**2*(-yy1 + yy2) + xx2**2*(yy1 - yy3) + xx1**2*(-yy2 + yy3))**2 &
544 : - 4*(xx3*(-yy1 + yy2) + xx2*(yy1 - yy3) + xx1*(-yy2 + yy3)) &
545 2 : *(xx1*xx3*(-xx1 + xx3)*yy2 + xx2**2*(xx3*yy1 - xx1*yy3) + xx2*(-xx3**2*yy1 + xx1**2*yy3))
546 2 : if (s2 < 0) then
547 0 : ierr = -1
548 0 : return
549 : end if
550 2 : b = sqrt(s2)
551 2 : a = xx3**2*(yy1 - yy2) + xx1**2*(yy2 - yy3) + xx2**2*(-yy1 + yy3)
552 2 : denom = 2*(xx3*(yy1 - yy2) + xx1*(yy2 - yy3) + xx2*(-yy1 + yy3))
553 2 : x = (a + b)/denom
554 2 : if (x > max(xx1, xx2, xx3)) x = (a - b)/denom
555 2 : if (x < min(xx1, xx2, xx3) .or. x > max(xx1, xx2, xx3)) ierr = -1
556 : end function find0_quadratic
557 :
558 1 : subroutine find_max_quadratic(x1, y1, x2, y2, x3, y3, xmax, ymax, ierr)
559 : ! x1 < x2 < x3 or x1 > x2 > x3. y2 > max(y1,y2)
560 : ! returns max location and value of quadratic fit to the three points
561 : real(dp), intent(in) :: x1, y1, x2, y2, x3, y3
562 : real(dp), intent(out) :: xmax, ymax
563 : integer, intent(out) :: ierr
564 1 : real(dp) :: a, b, c, dx1, dx2, dxmax
565 1 : ierr = 0; xmax = 0; ymax = 0
566 1 : dx1 = x2 - x1
567 1 : dx2 = x3 - x2
568 : ! f(dx) = a + b dx + c dx^2; dx = x - x1
569 1 : a = y1
570 1 : b = (y2 - y1)/dx1 + (y2 - y1)/(dx1 + dx2) + (dx1/dx2)*(y2 - y3)/(dx1 + dx2)
571 1 : c = (dx2*y1 - dx1*y2 - dx2*y2 + dx1*y3)/(dx1*dx2*(dx1 + dx2))
572 1 : dxmax = -b/(2d0*c)
573 1 : xmax = dxmax + x1
574 1 : ymax = a + dxmax*(b + dxmax*c)
575 1 : if (ymax < y2) ierr = -1
576 1 : if (y2 < max(y1, y2)) ierr = -1
577 1 : if (.not. ((x1 < x2 .and. x2 < x3) .or. (x1 > x2 .and. x2 > x3))) ierr = -1
578 1 : end subroutine find_max_quadratic
579 :
580 0 : subroutine two_piece_linear_coeffs(x, x0, x1, x2, a0, a1, a2, ierr)
581 : ! interpolation value at x is a0*f(x0) + a1*f(x1) + a2*f(x2)
582 : real(dp), intent(in) :: x, x0, x1, x2
583 : real(dp), intent(out) :: a0, a1, a2
584 : integer, intent(out) :: ierr
585 0 : ierr = 0
586 0 : if (x0 < x1 .and. x1 < x2) then
587 0 : if (x <= x0) then
588 0 : a0 = 1; a1 = 0; a2 = 0
589 0 : else if (x >= x2) then
590 0 : a0 = 0; a1 = 0; a2 = 1
591 0 : else if (x <= x1) then
592 0 : a1 = min(1d0, max(0d0, (x - x0)/(x1 - x0)))
593 0 : a0 = 1 - a1; a2 = 0
594 0 : else if (x < x2) then
595 0 : a2 = min(1d0, max(0d0, (x - x1)/(x2 - x1))) ! a2 => 1 as x => x2
596 0 : a1 = 1 - a2; a0 = 0
597 : end if
598 0 : else if (x0 > x1 .and. x1 > x2) then
599 0 : if (x >= x0) then
600 0 : a0 = 1; a1 = 0; a2 = 0
601 0 : else if (x <= x2) then
602 0 : a0 = 0; a1 = 0; a2 = 1
603 0 : else if (x >= x1) then
604 0 : a1 = min(1d0, max(0d0, (x - x0)/(x1 - x0)))
605 0 : a0 = 1 - a1; a2 = 0
606 0 : else if (x > x2) then
607 0 : a2 = min(1d0, max(0d0, (x - x1)/(x2 - x1))) ! a2 => 1 as x => x2
608 0 : a1 = 1 - a2; a0 = 0
609 : end if
610 : else
611 0 : ierr = -1
612 : end if
613 0 : end subroutine two_piece_linear_coeffs
614 :
615 7 : real(dp) function integrate(func, minx, maxx, args, atol, rtol, max_steps, ierr)
616 : procedure(integrand) :: func
617 : real(dp), intent(in) :: minx, maxx ! Min and max values to integrate over
618 : real(dp), intent(in) :: args(:) ! Extra args passed to func
619 : real(dp), intent(in) :: atol, rtol ! Absoulate and relative tolerances
620 : integer, intent(in) :: max_steps ! Max number of sub-steps
621 : integer, intent(inout) :: ierr ! Error code
622 7 : ierr = 0
623 7 : integrate = integrator(func, minx, maxx, args, atol, rtol, max_steps, ierr)
624 7 : end function integrate
625 :
626 : end module num_lib
|
