1# calculus.tcl -- 2# Package that implements several basic numerical methods, such 3# as the integration of a one-dimensional function and the 4# solution of a system of first-order differential equations. 5# 6# Copyright (c) 2002, 2003, 2004, 2006 by Arjen Markus. 7# Copyright (c) 2004 by Kevin B. Kenny. All rights reserved. 8# See the file "license.terms" for information on usage and redistribution 9# of this file, and for a DISCLAIMER OF ALL WARRANTIES. 10# 11# RCS: @(#) $Id: calculus.tcl,v 1.15 2008/10/08 03:30:48 andreas_kupries Exp $ 12 13package require Tcl 8.4 14package require math::interpolate 15package provide math::calculus 0.7.1 16 17# math::calculus -- 18# Namespace for the commands 19 20namespace eval ::math::calculus { 21 22 namespace import ::math::interpolate::neville 23 24 namespace import ::math::expectDouble ::math::expectInteger 25 26 namespace export \ 27 integral integralExpr integral2D integral3D \ 28 eulerStep heunStep rungeKuttaStep \ 29 boundaryValueSecondOrder solveTriDiagonal \ 30 newtonRaphson newtonRaphsonParameters 31 namespace export \ 32 integral2D_2accurate integral3D_accurate 33 34 namespace export romberg romberg_infinity 35 namespace export romberg_sqrtSingLower romberg_sqrtSingUpper 36 namespace export romberg_powerLawLower romberg_powerLawUpper 37 namespace export romberg_expLower romberg_expUpper 38 39 namespace export regula_falsi 40 41 variable nr_maxiter 20 42 variable nr_tolerance 0.001 43 44} 45 46# integral -- 47# Integrate a function over a given interval using the Simpson rule 48# 49# Arguments: 50# begin Start of the interval 51# end End of the interval 52# nosteps Number of steps in which to divide the interval 53# func Name of the function to be integrated (takes one 54# argument) 55# Return value: 56# Computed integral 57# 58proc ::math::calculus::integral { begin end nosteps func } { 59 60 set delta [expr {($end-$begin)/double($nosteps)}] 61 set hdelta [expr {$delta/2.0}] 62 set result 0.0 63 set xval $begin 64 set func_end [uplevel 1 $func $xval] 65 for { set i 1 } { $i <= $nosteps } { incr i } { 66 set func_begin $func_end 67 set func_middle [uplevel 1 $func [expr {$xval+$hdelta}]] 68 set func_end [uplevel 1 $func [expr {$xval+$delta}]] 69 set result [expr {$result+$func_begin+4.0*$func_middle+$func_end}] 70 71 set xval [expr {$begin+double($i)*$delta}] 72 } 73 74 return [expr {$result*$delta/6.0}] 75} 76 77# integralExpr -- 78# Integrate an expression with "x" as the integrate according to the 79# Simpson rule 80# 81# Arguments: 82# begin Start of the interval 83# end End of the interval 84# nosteps Number of steps in which to divide the interval 85# expression Expression with "x" as the integrate 86# Return value: 87# Computed integral 88# 89proc ::math::calculus::integralExpr { begin end nosteps expression } { 90 91 set delta [expr {($end-$begin)/double($nosteps)}] 92 set hdelta [expr {$delta/2.0}] 93 set result 0.0 94 set x $begin 95 # FRINK: nocheck 96 set func_end [expr $expression] 97 for { set i 1 } { $i <= $nosteps } { incr i } { 98 set func_begin $func_end 99 set x [expr {$x+$hdelta}] 100 # FRINK: nocheck 101 set func_middle [expr $expression] 102 set x [expr {$x+$hdelta}] 103 # FRINK: nocheck 104 set func_end [expr $expression] 105 set result [expr {$result+$func_begin+4.0*$func_middle+$func_end}] 106 107 set x [expr {$begin+double($i)*$delta}] 108 } 109 110 return [expr {$result*$delta/6.0}] 111} 112 113# integral2D -- 114# Integrate a given fucntion of two variables over a block, 115# using bilinear interpolation (for this moment: block function) 116# 117# Arguments: 118# xinterval Start, stop and number of steps of the "x" interval 119# yinterval Start, stop and number of steps of the "y" interval 120# func Function of the two variables to be integrated 121# Return value: 122# Computed integral 123# 124proc ::math::calculus::integral2D { xinterval yinterval func } { 125 126 foreach { xbegin xend xnumber } $xinterval { break } 127 foreach { ybegin yend ynumber } $yinterval { break } 128 129 set xdelta [expr {($xend-$xbegin)/double($xnumber)}] 130 set ydelta [expr {($yend-$ybegin)/double($ynumber)}] 131 set hxdelta [expr {$xdelta/2.0}] 132 set hydelta [expr {$ydelta/2.0}] 133 set result 0.0 134 set dxdy [expr {$xdelta*$ydelta}] 135 for { set j 0 } { $j < $ynumber } { incr j } { 136 set y [expr {$ybegin+$hydelta+double($j)*$ydelta}] 137 for { set i 0 } { $i < $xnumber } { incr i } { 138 set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}] 139 set func_value [uplevel 1 $func $x $y] 140 set result [expr {$result+$func_value}] 141 } 142 } 143 144 return [expr {$result*$dxdy}] 145} 146 147# integral3D -- 148# Integrate a given fucntion of two variables over a block, 149# using trilinear interpolation (for this moment: block function) 150# 151# Arguments: 152# xinterval Start, stop and number of steps of the "x" interval 153# yinterval Start, stop and number of steps of the "y" interval 154# zinterval Start, stop and number of steps of the "z" interval 155# func Function of the three variables to be integrated 156# Return value: 157# Computed integral 158# 159proc ::math::calculus::integral3D { xinterval yinterval zinterval func } { 160 161 foreach { xbegin xend xnumber } $xinterval { break } 162 foreach { ybegin yend ynumber } $yinterval { break } 163 foreach { zbegin zend znumber } $zinterval { break } 164 165 set xdelta [expr {($xend-$xbegin)/double($xnumber)}] 166 set ydelta [expr {($yend-$ybegin)/double($ynumber)}] 167 set zdelta [expr {($zend-$zbegin)/double($znumber)}] 168 set hxdelta [expr {$xdelta/2.0}] 169 set hydelta [expr {$ydelta/2.0}] 170 set hzdelta [expr {$zdelta/2.0}] 171 set result 0.0 172 set dxdydz [expr {$xdelta*$ydelta*$zdelta}] 173 for { set k 0 } { $k < $znumber } { incr k } { 174 set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}] 175 for { set j 0 } { $j < $ynumber } { incr j } { 176 set y [expr {$ybegin+$hydelta+double($j)*$ydelta}] 177 for { set i 0 } { $i < $xnumber } { incr i } { 178 set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}] 179 set func_value [uplevel 1 $func $x $y $z] 180 set result [expr {$result+$func_value}] 181 } 182 } 183 } 184 185 return [expr {$result*$dxdydz}] 186} 187 188# integral2D_accurate -- 189# Integrate a given function of two variables over a block, 190# using a four-point quadrature formula 191# 192# Arguments: 193# xinterval Start, stop and number of steps of the "x" interval 194# yinterval Start, stop and number of steps of the "y" interval 195# func Function of the two variables to be integrated 196# Return value: 197# Computed integral 198# 199proc ::math::calculus::integral2D_accurate { xinterval yinterval func } { 200 201 foreach { xbegin xend xnumber } $xinterval { break } 202 foreach { ybegin yend ynumber } $yinterval { break } 203 204 set alpha [expr {sqrt(2.0/3.0)}] 205 set minalpha [expr {-$alpha}] 206 set dpoints [list $alpha 0.0 $minalpha 0.0 0.0 $alpha 0.0 $minalpha] 207 208 set xdelta [expr {($xend-$xbegin)/double($xnumber)}] 209 set ydelta [expr {($yend-$ybegin)/double($ynumber)}] 210 set hxdelta [expr {$xdelta/2.0}] 211 set hydelta [expr {$ydelta/2.0}] 212 set result 0.0 213 set dxdy [expr {0.25*$xdelta*$ydelta}] 214 215 for { set j 0 } { $j < $ynumber } { incr j } { 216 set y [expr {$ybegin+$hydelta+double($j)*$ydelta}] 217 for { set i 0 } { $i < $xnumber } { incr i } { 218 set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}] 219 220 foreach {dx dy} $dpoints { 221 set x1 [expr {$x+$dx}] 222 set y1 [expr {$y+$dy}] 223 set func_value [uplevel 1 $func $x1 $y1] 224 set result [expr {$result+$func_value}] 225 } 226 } 227 } 228 229 return [expr {$result*$dxdy}] 230} 231 232# integral3D_accurate -- 233# Integrate a given function of three variables over a block, 234# using an 8-point quadrature formula 235# 236# Arguments: 237# xinterval Start, stop and number of steps of the "x" interval 238# yinterval Start, stop and number of steps of the "y" interval 239# zinterval Start, stop and number of steps of the "z" interval 240# func Function of the three variables to be integrated 241# Return value: 242# Computed integral 243# 244proc ::math::calculus::integral3D_accurate { xinterval yinterval zinterval func } { 245 246 foreach { xbegin xend xnumber } $xinterval { break } 247 foreach { ybegin yend ynumber } $yinterval { break } 248 foreach { zbegin zend znumber } $zinterval { break } 249 250 set alpha [expr {sqrt(1.0/3.0)}] 251 set minalpha [expr {-$alpha}] 252 253 set dpoints [list $alpha $alpha $alpha \ 254 $alpha $alpha $minalpha \ 255 $alpha $minalpha $alpha \ 256 $alpha $minalpha $minalpha \ 257 $minalpha $alpha $alpha \ 258 $minalpha $alpha $minalpha \ 259 $minalpha $minalpha $alpha \ 260 $minalpha $minalpha $minalpha ] 261 262 set xdelta [expr {($xend-$xbegin)/double($xnumber)}] 263 set ydelta [expr {($yend-$ybegin)/double($ynumber)}] 264 set zdelta [expr {($zend-$zbegin)/double($znumber)}] 265 set hxdelta [expr {$xdelta/2.0}] 266 set hydelta [expr {$ydelta/2.0}] 267 set hzdelta [expr {$zdelta/2.0}] 268 set result 0.0 269 set dxdydz [expr {0.125*$xdelta*$ydelta*$zdelta}] 270 271 for { set k 0 } { $k < $znumber } { incr k } { 272 set z [expr {$zbegin+$hzdelta+double($k)*$zdelta}] 273 for { set j 0 } { $j < $ynumber } { incr j } { 274 set y [expr {$ybegin+$hydelta+double($j)*$ydelta}] 275 for { set i 0 } { $i < $xnumber } { incr i } { 276 set x [expr {$xbegin+$hxdelta+double($i)*$xdelta}] 277 278 foreach {dx dy dz} $dpoints { 279 set x1 [expr {$x+$dx}] 280 set y1 [expr {$y+$dy}] 281 set z1 [expr {$z+$dz}] 282 set func_value [uplevel 1 $func $x1 $y1 $z1] 283 set result [expr {$result+$func_value}] 284 } 285 } 286 } 287 } 288 289 return [expr {$result*$dxdydz}] 290} 291 292# eulerStep -- 293# Integrate a system of ordinary differential equations of the type 294# x' = f(x,t), where x is a vector of quantities. Integration is 295# done over a single step according to Euler's method. 296# 297# Arguments: 298# t Start value of independent variable (time for instance) 299# tstep Step size of interval 300# xvec Vector of dependent values at the start 301# func Function taking the arguments t and xvec to return 302# the derivative of each dependent variable. 303# Return value: 304# List of values at the end of the step 305# 306proc ::math::calculus::eulerStep { t tstep xvec func } { 307 308 set xderiv [uplevel 1 $func $t [list $xvec]] 309 set result {} 310 foreach xv $xvec dx $xderiv { 311 set xnew [expr {$xv+$tstep*$dx}] 312 lappend result $xnew 313 } 314 315 return $result 316} 317 318# heunStep -- 319# Integrate a system of ordinary differential equations of the type 320# x' = f(x,t), where x is a vector of quantities. Integration is 321# done over a single step according to Heun's method. 322# 323# Arguments: 324# t Start value of independent variable (time for instance) 325# tstep Step size of interval 326# xvec Vector of dependent values at the start 327# func Function taking the arguments t and xvec to return 328# the derivative of each dependent variable. 329# Return value: 330# List of values at the end of the step 331# 332proc ::math::calculus::heunStep { t tstep xvec func } { 333 334 # 335 # Predictor step 336 # 337 set funcq [uplevel 1 namespace which -command $func] 338 set xpred [eulerStep $t $tstep $xvec $funcq] 339 340 # 341 # Corrector step 342 # 343 set tcorr [expr {$t+$tstep}] 344 set xcorr [eulerStep $t $tstep $xpred $funcq] 345 346 set result {} 347 foreach xv $xvec xc $xcorr { 348 set xnew [expr {0.5*($xv+$xc)}] 349 lappend result $xnew 350 } 351 352 return $result 353} 354 355# rungeKuttaStep -- 356# Integrate a system of ordinary differential equations of the type 357# x' = f(x,t), where x is a vector of quantities. Integration is 358# done over a single step according to Runge-Kutta 4th order. 359# 360# Arguments: 361# t Start value of independent variable (time for instance) 362# tstep Step size of interval 363# xvec Vector of dependent values at the start 364# func Function taking the arguments t and xvec to return 365# the derivative of each dependent variable. 366# Return value: 367# List of values at the end of the step 368# 369proc ::math::calculus::rungeKuttaStep { t tstep xvec func } { 370 371 set funcq [uplevel 1 namespace which -command $func] 372 373 # 374 # Four steps: 375 # - k1 = tstep*func(t,x0) 376 # - k2 = tstep*func(t+0.5*tstep,x0+0.5*k1) 377 # - k3 = tstep*func(t+0.5*tstep,x0+0.5*k2) 378 # - k4 = tstep*func(t+ tstep,x0+ k3) 379 # - x1 = x0 + (k1+2*k2+2*k3+k4)/6 380 # 381 set tstep2 [expr {$tstep/2.0}] 382 set tstep6 [expr {$tstep/6.0}] 383 384 set xk1 [$funcq $t $xvec] 385 set xvec2 {} 386 foreach x1 $xvec xv $xk1 { 387 lappend xvec2 [expr {$x1+$tstep2*$xv}] 388 } 389 set xk2 [$funcq [expr {$t+$tstep2}] $xvec2] 390 391 set xvec3 {} 392 foreach x1 $xvec xv $xk2 { 393 lappend xvec3 [expr {$x1+$tstep2*$xv}] 394 } 395 set xk3 [$funcq [expr {$t+$tstep2}] $xvec3] 396 397 set xvec4 {} 398 foreach x1 $xvec xv $xk3 { 399 lappend xvec4 [expr {$x1+$tstep*$xv}] 400 } 401 set xk4 [$funcq [expr {$t+$tstep}] $xvec4] 402 403 set result {} 404 foreach x0 $xvec k1 $xk1 k2 $xk2 k3 $xk3 k4 $xk4 { 405 set dx [expr {$k1+2.0*$k2+2.0*$k3+$k4}] 406 lappend result [expr {$x0+$dx*$tstep6}] 407 } 408 409 return $result 410} 411 412# boundaryValueSecondOrder -- 413# Integrate a second-order differential equation and solve for 414# given boundary values. 415# 416# The equation is (see the documentation): 417# d dy d 418# -- A(x) -- + -- B(x) y + C(x) y = D(x) 419# dx dx dx 420# 421# The procedure uses finite differences and tridiagonal matrices to 422# solve the equation. The boundary values are put in the matrix 423# directly. 424# 425# Arguments: 426# coeff_func Name of triple-valued function for coefficients A, B, C 427# force_func Name of the function providing the force term D(x) 428# leftbnd Left boundary condition (list of: xvalue, boundary 429# value or keyword zero-flux, zero-derivative) 430# rightbnd Right boundary condition (ditto) 431# nostep Number of steps 432# Return value: 433# List of x-values and calculated values (x1, y1, x2, y2, ...) 434# 435proc ::math::calculus::boundaryValueSecondOrder { 436 coeff_func force_func leftbnd rightbnd nostep } { 437 438 set coeffq [uplevel 1 namespace which -command $coeff_func] 439 set forceq [uplevel 1 namespace which -command $force_func] 440 441 if { [llength $leftbnd] != 2 || [llength $rightbnd] != 2 } { 442 error "Boundary condition(s) incorrect" 443 } 444 if { $nostep < 1 } { 445 error "Number of steps must be larger/equal 1" 446 } 447 448 # 449 # Set up the matrix, as three different lists and the 450 # righthand side as the fourth 451 # 452 set xleft [lindex $leftbnd 0] 453 set xright [lindex $rightbnd 0] 454 set xstep [expr {($xright-$xleft)/double($nostep)}] 455 456 set acoeff {} 457 set bcoeff {} 458 set ccoeff {} 459 set dvalue {} 460 461 set x $xleft 462 foreach {A B C} [$coeffq $x] { break } 463 464 set A1 [expr {$A/$xstep-0.5*$B}] 465 set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}] 466 set C1 0.0 467 468 for { set i 1 } { $i <= $nostep } { incr i } { 469 set x [expr {$xleft+double($i)*$xstep}] 470 if { [expr {abs($x)-0.5*abs($xstep)}] < 0.0 } { 471 set x 0.0 472 } 473 foreach {A B C} [$coeffq $x] { break } 474 475 set A2 0.0 476 set B2 [expr {$A/$xstep-0.5*$B+0.5*$C*$xstep}] 477 set C2 [expr {$A/$xstep+0.5*$B}] 478 lappend acoeff [expr {$A1+$A2}] 479 lappend bcoeff [expr {-$B1-$B2}] 480 lappend ccoeff [expr {$C1+$C2}] 481 set A1 [expr {$A/$xstep-0.5*$B}] 482 set B1 [expr {$A/$xstep+0.5*$B+0.5*$C*$xstep}] 483 set C1 0.0 484 } 485 set xvec {} 486 for { set i 0 } { $i < $nostep } { incr i } { 487 set x [expr {$xleft+(0.5+double($i))*$xstep}] 488 if { [expr {abs($x)-0.25*abs($xstep)}] < 0.0 } { 489 set x 0.0 490 } 491 lappend xvec $x 492 lappend dvalue [expr {$xstep*[$forceq $x]}] 493 } 494 495 # 496 # Substitute the boundary values 497 # 498 set A [lindex $acoeff 0] 499 set D [lindex $dvalue 0] 500 set D1 [expr {$D-$A*[lindex $leftbnd 1]}] 501 set C [lindex $ccoeff end] 502 set D [lindex $dvalue end] 503 set D2 [expr {$D-$C*[lindex $rightbnd 1]}] 504 set dvalue [concat $D1 [lrange $dvalue 1 end-1] $D2] 505 506 set yvec [solveTriDiagonal [lrange $acoeff 1 end] $bcoeff [lrange $ccoeff 0 end-1] $dvalue] 507 508 foreach x $xvec y $yvec { 509 lappend result $x $y 510 } 511 return $result 512} 513 514# solveTriDiagonal -- 515# Solve a system of equations Ax = b where A is a tridiagonal matrix 516# 517# Arguments: 518# acoeff Values on lower diagonal 519# bcoeff Values on main diagonal 520# ccoeff Values on upper diagonal 521# dvalue Values on righthand side 522# Return value: 523# List of values forming the solution 524# 525proc ::math::calculus::solveTriDiagonal { acoeff bcoeff ccoeff dvalue } { 526 527 set nostep [llength $acoeff] 528 # 529 # First step: Gauss-elimination 530 # 531 set B [lindex $bcoeff 0] 532 set C [lindex $ccoeff 0] 533 set D [lindex $dvalue 0] 534 set acoeff [concat 0.0 $acoeff] 535 set bcoeff2 [list $B] 536 set dvalue2 [list $D] 537 for { set i 1 } { $i <= $nostep } { incr i } { 538 set A2 [lindex $acoeff $i] 539 set B2 [lindex $bcoeff $i] 540 set D2 [lindex $dvalue $i] 541 set ratab [expr {$A2/double($B)}] 542 set B2 [expr {$B2-$ratab*$C}] 543 set D2 [expr {$D2-$ratab*$D}] 544 lappend bcoeff2 $B2 545 lappend dvalue2 $D2 546 set B $B2 547 set C [lindex $ccoeff $i] 548 set D $D2 549 } 550 551 # 552 # Second step: substitution 553 # 554 set yvec {} 555 set B [lindex $bcoeff2 end] 556 set D [lindex $dvalue2 end] 557 set y [expr {$D/$B}] 558 for { set i [expr {$nostep-1}] } { $i >= 0 } { incr i -1 } { 559 set yvec [concat $y $yvec] 560 set B [lindex $bcoeff2 $i] 561 set C [lindex $ccoeff $i] 562 set D [lindex $dvalue2 $i] 563 set y [expr {($D-$C*$y)/$B}] 564 } 565 set yvec [concat $y $yvec] 566 567 return $yvec 568} 569 570# newtonRaphson -- 571# Determine the root of an equation via the Newton-Raphson method 572# 573# Arguments: 574# func Function (proc) in x 575# deriv Derivative (proc) of func w.r.t. x 576# initval Initial value for x 577# Return value: 578# Estimate of root 579# 580proc ::math::calculus::newtonRaphson { func deriv initval } { 581 variable nr_maxiter 582 variable nr_tolerance 583 584 set funcq [uplevel 1 namespace which -command $func] 585 set derivq [uplevel 1 namespace which -command $deriv] 586 587 set value $initval 588 set diff [expr {10.0*$nr_tolerance}] 589 590 for { set i 0 } { $i < $nr_maxiter } { incr i } { 591 if { $diff < $nr_tolerance } { 592 break 593 } 594 595 set newval [expr {$value-[$funcq $value]/[$derivq $value]}] 596 if { $value != 0.0 } { 597 set diff [expr {abs($newval-$value)/abs($value)}] 598 } else { 599 set diff [expr {abs($newval-$value)}] 600 } 601 set value $newval 602 } 603 604 return $newval 605} 606 607# newtonRaphsonParameters -- 608# Set the parameters for the Newton-Raphson method 609# 610# Arguments: 611# maxiter Maximum number of iterations 612# tolerance Relative precisiion of the result 613# Return value: 614# None 615# 616proc ::math::calculus::newtonRaphsonParameters { maxiter tolerance } { 617 variable nr_maxiter 618 variable nr_tolerance 619 620 if { $maxiter > 0 } { 621 set nr_maxiter $maxiter 622 } 623 if { $tolerance > 0 } { 624 set nr_tolerance $tolerance 625 } 626} 627 628#---------------------------------------------------------------------- 629# 630# midpoint -- 631# 632# Perform one set of steps in evaluating an integral using the 633# midpoint method. 634# 635# Usage: 636# midpoint f a b s ?n? 637# 638# Parameters: 639# f - function to integrate 640# a - One limit of integration 641# b - Other limit of integration. a and b need not be in ascending 642# order. 643# s - Value returned from a previous call to midpoint (see below) 644# n - Step number (see below) 645# 646# Results: 647# Returns an estimate of the integral obtained by dividing the 648# interval into 3**n equal intervals and using the midpoint rule. 649# 650# Side effects: 651# f is evaluated 2*3**(n-1) times and may have side effects. 652# 653# The 'midpoint' procedure is designed for successive approximations. 654# It should be called initially with n==0. On this initial call, s 655# is ignored. The function is evaluated at the midpoint of the interval, and 656# the value is multiplied by the width of the interval to give the 657# coarsest possible estimate of the integral. 658# 659# On each iteration except the first, n should be incremented by one, 660# and the previous value returned from [midpoint] should be supplied 661# as 's'. The function will be evaluated at additional points 662# to give a total of 3**n equally spaced points, and the estimate 663# of the integral will be updated and returned 664# 665# Under normal circumstances, user code will not call this function 666# directly. Instead, it will use ::math::calculus::romberg to 667# do error control and extrapolation to a zero step size. 668# 669#---------------------------------------------------------------------- 670 671proc ::math::calculus::midpoint { f a b { n 0 } { s 0. } } { 672 673 if { $n == 0 } { 674 675 # First iteration. Simply evaluate the function at the midpoint 676 # of the interval. 677 678 set cmd $f; lappend cmd [expr { 0.5 * ( $a + $b ) }]; set v [eval $cmd] 679 return [expr { ( $b - $a ) * $v }] 680 681 } else { 682 683 # Subsequent iterations. We've divided the interval into 684 # $it subintervals. Evaluate the function at the 1/3 and 685 # 2/3 points of each subinterval. Then update the estimate 686 # of the integral that we produced on the last step with 687 # the new sum. 688 689 set it [expr { pow( 3, $n-1 ) }] 690 set h [expr { ( $b - $a ) / ( 3. * $it ) }] 691 set h2 [expr { $h + $h }] 692 set x [expr { $a + 0.5 * $h }] 693 set sum 0 694 for { set j 0 } { $j < $it } { incr j } { 695 set cmd $f; lappend cmd $x; set y [eval $cmd] 696 set sum [expr { $sum + $y }] 697 set x [expr { $x + $h2 }] 698 set cmd $f; lappend cmd $x; set y [eval $cmd] 699 set sum [expr { $sum + $y }] 700 set x [expr { $x + $h}] 701 } 702 return [expr { ( $s + ( $b - $a ) * $sum / $it ) / 3. }] 703 704 } 705} 706 707#---------------------------------------------------------------------- 708# 709# romberg -- 710# 711# Compute the integral of a function over an interval using 712# Romberg's method. 713# 714# Usage: 715# romberg f a b ?-option value?... 716# 717# Parameters: 718# f - Function to integrate. Must be a single Tcl command, 719# to which will be appended the abscissa at which the function 720# should be evaluated. f should be analytic over the 721# region of integration, but may have a removable singularity 722# at either endpoint. 723# a - One bound of the interval 724# b - The other bound of the interval. a and b need not be in 725# ascending order. 726# 727# Options: 728# -abserror ABSERROR 729# Requests that the integration be performed to make 730# the estimated absolute error of the integral less than 731# the given value. Default is 1.e-10. 732# -relerror RELERROR 733# Requests that the integration be performed to make 734# the estimated absolute error of the integral less than 735# the given value. Default is 1.e-6. 736# -degree N 737# Specifies the degree of the polynomial that will be 738# used to extrapolate to a zero step size. -degree 0 739# requests integration with the midpoint rule; -degree 1 740# is equivalent to Simpson's 3/8 rule; higher degrees 741# are difficult to describe but (within reason) give 742# faster convergence for smooth functions. Default is 743# -degree 4. 744# -maxiter N 745# Specifies the maximum number of triplings of the 746# number of steps to take in integration. At most 747# 3**N function evaluations will be performed in 748# integrating with -maxiter N. The integration 749# will terminate at that time, even if the result 750# satisfies neither the -relerror nor -abserror tests. 751# 752# Results: 753# Returns a two-element list. The first element is the estimated 754# value of the integral; the second is the estimated absolute 755# error of the value. 756# 757#---------------------------------------------------------------------- 758 759proc ::math::calculus::romberg { f a b args } { 760 761 # Replace f with a context-independent version 762 763 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 764 765 # Assign default parameters 766 767 array set params { 768 -abserror 1.0e-10 769 -degree 4 770 -relerror 1.0e-6 771 -maxiter 14 772 } 773 774 # Extract parameters 775 776 if { ( [llength $args] % 2 ) != 0 } { 777 return -code error -errorcode [list romberg wrongNumArgs] \ 778 "wrong \# args, should be\ 779 \"[lreplace [info level 0] 1 end \ 780 f x1 x2 ?-option value?...]\"" 781 } 782 foreach { key value } $args { 783 if { ![info exists params($key)] } { 784 return -code error -errorcode [list romberg badoption $key] \ 785 "unknown option \"$key\",\ 786 should be -abserror, -degree, -relerror, or -maxiter" 787 } 788 set params($key) $value 789 } 790 791 # Check params 792 793 if { ![string is double -strict $a] } { 794 return -code error [expectDouble $a] 795 } 796 if { ![string is double -strict $b] } { 797 return -code error [expectDouble $b] 798 } 799 if { ![string is double -strict $params(-abserror)] } { 800 return -code error [expectDouble $params(-abserror)] 801 } 802 if { ![string is integer -strict $params(-degree)] } { 803 return -code error [expectInteger $params(-degree)] 804 } 805 if { ![string is integer -strict $params(-maxiter)] } { 806 return -code error [expectInteger $params(-maxiter)] 807 } 808 if { ![string is double -strict $params(-relerror)] } { 809 return -code error [expectDouble $params(-relerror)] 810 } 811 foreach key {-abserror -degree -maxiter -relerror} { 812 if { $params($key) <= 0 } { 813 return -code error -errorcode [list romberg notPositive $key] \ 814 "$key must be positive" 815 } 816 } 817 if { $params(-maxiter) <= $params(-degree) } { 818 return -code error -errorcode [list romberg tooFewIter] \ 819 "-maxiter must be greater than -degree" 820 } 821 822 # Create lists of step size and sum with the given number of steps. 823 824 set x [list] 825 set y [list] 826 set s 0; # Current best estimate of integral 827 set indx end-$params(-degree) 828 set pow3 1.; # Current step size (times b-a) 829 830 # Perform successive integrations, tripling the number of steps each time 831 832 for { set i 0 } { $i < $params(-maxiter) } { incr i } { 833 set s [midpoint $f $a $b $i $s] 834 lappend x $pow3 835 lappend y $s 836 set pow3 [expr { $pow3 / 9. }] 837 838 # Once $degree steps have been done, start Richardson extrapolation 839 # to a zero step size. 840 841 if { $i >= $params(-degree) } { 842 set x [lrange $x $indx end] 843 set y [lrange $y $indx end] 844 foreach {estimate err} [neville $x $y 0.] break 845 if { $err < $params(-abserror) 846 || $err < $params(-relerror) * abs($estimate) } { 847 return [list $estimate $err] 848 } 849 } 850 } 851 852 # If -maxiter iterations have been done, give up, and return 853 # with the current error estimate. 854 855 return [list $estimate $err] 856} 857 858#---------------------------------------------------------------------- 859# 860# u_infinity -- 861# Change of variable for integrating over a half-infinite 862# interval 863# 864# Parameters: 865# f - Function being integrated 866# u - 1/x, where x is the abscissa where f is to be evaluated 867# 868# Results: 869# Returns f(1/u)/(u**2) 870# 871# Side effects: 872# Whatever f does. 873# 874#---------------------------------------------------------------------- 875 876proc ::math::calculus::u_infinity { f u } { 877 set cmd $f 878 lappend cmd [expr { 1.0 / $u }] 879 set y [eval $cmd] 880 return [expr { $y / ( $u * $u ) }] 881} 882 883#---------------------------------------------------------------------- 884# 885# romberg_infinity -- 886# Evaluate a function on a half-open interval 887# 888# Usage: 889# Same as 'romberg' 890# 891# The 'romberg_infinity' procedure performs Romberg integration on 892# an interval [a,b] where an infinite a or b may be represented by 893# a large number (e.g. 1.e30). It operates by a change of variable; 894# instead of integrating f(x) from a to b, it makes a change 895# of variable u = 1/x, and integrates from 1/b to 1/a f(1/u)/u**2 du. 896# 897#---------------------------------------------------------------------- 898 899proc ::math::calculus::romberg_infinity { f a b args } { 900 if { ![string is double -strict $a] } { 901 return -code error [expectDouble $a] 902 } 903 if { ![string is double -strict $b] } { 904 return -code error [expectDouble $b] 905 } 906 if { $a * $b <= 0. } { 907 return -code error -errorcode {romberg_infinity cross-axis} \ 908 "limits of integration have opposite sign" 909 } 910 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 911 set f [list u_infinity $f] 912 return [eval [linsert $args 0 \ 913 romberg $f [expr { 1.0 / $b }] [expr { 1.0 / $a }]]] 914} 915 916#---------------------------------------------------------------------- 917# 918# u_sqrtSingLower -- 919# Change of variable for integrating over an interval with 920# an inverse square root singularity at the lower bound. 921# 922# Parameters: 923# f - Function being integrated 924# a - Lower bound 925# u - sqrt(x-a), where x is the abscissa where f is to be evaluated 926# 927# Results: 928# Returns 2 * u * f( a + u**2 ) 929# 930# Side effects: 931# Whatever f does. 932# 933#---------------------------------------------------------------------- 934 935proc ::math::calculus::u_sqrtSingLower { f a u } { 936 set cmd $f 937 lappend cmd [expr { $a + $u * $u }] 938 set y [eval $cmd] 939 return [expr { 2. * $u * $y }] 940} 941 942#---------------------------------------------------------------------- 943# 944# u_sqrtSingUpper -- 945# Change of variable for integrating over an interval with 946# an inverse square root singularity at the upper bound. 947# 948# Parameters: 949# f - Function being integrated 950# b - Upper bound 951# u - sqrt(b-x), where x is the abscissa where f is to be evaluated 952# 953# Results: 954# Returns 2 * u * f( b - u**2 ) 955# 956# Side effects: 957# Whatever f does. 958# 959#---------------------------------------------------------------------- 960 961proc ::math::calculus::u_sqrtSingUpper { f b u } { 962 set cmd $f 963 lappend cmd [expr { $b - $u * $u }] 964 set y [eval $cmd] 965 return [expr { 2. * $u * $y }] 966} 967 968#---------------------------------------------------------------------- 969# 970# math::calculus::romberg_sqrtSingLower -- 971# Integrate a function with an inverse square root singularity 972# at the lower bound 973# 974# Usage: 975# Same as 'romberg' 976# 977# The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg' 978# for integrating a function with an inverse square root singularity 979# at the lower bound of the interval. It works by making the change 980# of variable u = sqrt( x-a ). 981# 982#---------------------------------------------------------------------- 983 984proc ::math::calculus::romberg_sqrtSingLower { f a b args } { 985 if { ![string is double -strict $a] } { 986 return -code error [expectDouble $a] 987 } 988 if { ![string is double -strict $b] } { 989 return -code error [expectDouble $b] 990 } 991 if { $a >= $b } { 992 return -code error "limits of integration out of order" 993 } 994 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 995 set f [list u_sqrtSingLower $f $a] 996 return [eval [linsert $args 0 \ 997 romberg $f 0 [expr { sqrt( $b - $a ) }]]] 998} 999 1000#---------------------------------------------------------------------- 1001# 1002# math::calculus::romberg_sqrtSingUpper -- 1003# Integrate a function with an inverse square root singularity 1004# at the upper bound 1005# 1006# Usage: 1007# Same as 'romberg' 1008# 1009# The 'romberg_sqrtSingUpper' procedure is a wrapper for 'romberg' 1010# for integrating a function with an inverse square root singularity 1011# at the upper bound of the interval. It works by making the change 1012# of variable u = sqrt( b-x ). 1013# 1014#---------------------------------------------------------------------- 1015 1016proc ::math::calculus::romberg_sqrtSingUpper { f a b args } { 1017 if { ![string is double -strict $a] } { 1018 return -code error [expectDouble $a] 1019 } 1020 if { ![string is double -strict $b] } { 1021 return -code error [expectDouble $b] 1022 } 1023 if { $a >= $b } { 1024 return -code error "limits of integration out of order" 1025 } 1026 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 1027 set f [list u_sqrtSingUpper $f $b] 1028 return [eval [linsert $args 0 \ 1029 romberg $f 0. [expr { sqrt( $b - $a ) }]]] 1030} 1031 1032#---------------------------------------------------------------------- 1033# 1034# u_powerLawLower -- 1035# Change of variable for integrating over an interval with 1036# an integrable power law singularity at the lower bound. 1037# 1038# Parameters: 1039# f - Function being integrated 1040# gammaover1mgamma - gamma / (1 - gamma), where gamma is the power 1041# oneover1mgamma - 1 / (1 - gamma), where gamma is the power 1042# a - Lower limit of integration 1043# u - Changed variable u == (x-a)**(1-gamma) 1044# 1045# Results: 1046# Returns u**(1/1-gamma) * f(a + u**(1/1-gamma) ). 1047# 1048# Side effects: 1049# Whatever f does. 1050# 1051#---------------------------------------------------------------------- 1052 1053proc ::math::calculus::u_powerLawLower { f gammaover1mgamma oneover1mgamma 1054 a u } { 1055 set cmd $f 1056 lappend cmd [expr { $a + pow( $u, $oneover1mgamma ) }] 1057 set y [eval $cmd] 1058 return [expr { $y * pow( $u, $gammaover1mgamma ) }] 1059} 1060 1061#---------------------------------------------------------------------- 1062# 1063# math::calculus::romberg_powerLawLower -- 1064# Integrate a function with an integrable power law singularity 1065# at the lower bound 1066# 1067# Usage: 1068# romberg_powerLawLower gamma f a b ?-option value...? 1069# 1070# Parameters: 1071# gamma - Power (0<gamma<1) of the singularity 1072# f - Function to integrate. Must be a single Tcl command, 1073# to which will be appended the abscissa at which the function 1074# should be evaluated. f is expected to have an integrable 1075# power law singularity at the lower endpoint; that is, the 1076# integrand is expected to diverge as (x-a)**gamma. 1077# a - One bound of the interval 1078# b - The other bound of the interval. a and b need not be in 1079# ascending order. 1080# 1081# Options: 1082# -abserror ABSERROR 1083# Requests that the integration be performed to make 1084# the estimated absolute error of the integral less than 1085# the given value. Default is 1.e-10. 1086# -relerror RELERROR 1087# Requests that the integration be performed to make 1088# the estimated absolute error of the integral less than 1089# the given value. Default is 1.e-6. 1090# -degree N 1091# Specifies the degree of the polynomial that will be 1092# used to extrapolate to a zero step size. -degree 0 1093# requests integration with the midpoint rule; -degree 1 1094# is equivalent to Simpson's 3/8 rule; higher degrees 1095# are difficult to describe but (within reason) give 1096# faster convergence for smooth functions. Default is 1097# -degree 4. 1098# -maxiter N 1099# Specifies the maximum number of triplings of the 1100# number of steps to take in integration. At most 1101# 3**N function evaluations will be performed in 1102# integrating with -maxiter N. The integration 1103# will terminate at that time, even if the result 1104# satisfies neither the -relerror nor -abserror tests. 1105# 1106# Results: 1107# Returns a two-element list. The first element is the estimated 1108# value of the integral; the second is the estimated absolute 1109# error of the value. 1110# 1111# The 'romberg_sqrtSingLower' procedure is a wrapper for 'romberg' 1112# for integrating a function with an integrable power law singularity 1113# at the lower bound of the interval. It works by making the change 1114# of variable u = (x-a)**(1-gamma). 1115# 1116#---------------------------------------------------------------------- 1117 1118proc ::math::calculus::romberg_powerLawLower { gamma f a b args } { 1119 if { ![string is double -strict $gamma] } { 1120 return -code error [expectDouble $gamma] 1121 } 1122 if { $gamma <= 0.0 || $gamma >= 1.0 } { 1123 return -code error -errorcode [list romberg gammaTooBig] \ 1124 "gamma must lie in the interval (0,1)" 1125 } 1126 if { ![string is double -strict $a] } { 1127 return -code error [expectDouble $a] 1128 } 1129 if { ![string is double -strict $b] } { 1130 return -code error [expectDouble $b] 1131 } 1132 if { $a >= $b } { 1133 return -code error "limits of integration out of order" 1134 } 1135 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 1136 set onemgamma [expr { 1. - $gamma }] 1137 set f [list u_powerLawLower $f \ 1138 [expr { $gamma / $onemgamma }] \ 1139 [expr { 1 / $onemgamma }] \ 1140 $a] 1141 1142 set limit [expr { pow( $b - $a, $onemgamma ) }] 1143 set result {} 1144 foreach v [eval [linsert $args 0 romberg $f 0 $limit]] { 1145 lappend result [expr { $v / $onemgamma }] 1146 } 1147 return $result 1148 1149} 1150 1151#---------------------------------------------------------------------- 1152# 1153# u_powerLawLower -- 1154# Change of variable for integrating over an interval with 1155# an integrable power law singularity at the upper bound. 1156# 1157# Parameters: 1158# f - Function being integrated 1159# gammaover1mgamma - gamma / (1 - gamma), where gamma is the power 1160# oneover1mgamma - 1 / (1 - gamma), where gamma is the power 1161# b - Upper limit of integration 1162# u - Changed variable u == (b-x)**(1-gamma) 1163# 1164# Results: 1165# Returns u**(1/1-gamma) * f(b-u**(1/1-gamma) ). 1166# 1167# Side effects: 1168# Whatever f does. 1169# 1170#---------------------------------------------------------------------- 1171 1172proc ::math::calculus::u_powerLawUpper { f gammaover1mgamma oneover1mgamma 1173 b u } { 1174 set cmd $f 1175 lappend cmd [expr { $b - pow( $u, $oneover1mgamma ) }] 1176 set y [eval $cmd] 1177 return [expr { $y * pow( $u, $gammaover1mgamma ) }] 1178} 1179 1180#---------------------------------------------------------------------- 1181# 1182# math::calculus::romberg_powerLawUpper -- 1183# Integrate a function with an integrable power law singularity 1184# at the upper bound 1185# 1186# Usage: 1187# romberg_powerLawLower gamma f a b ?-option value...? 1188# 1189# Parameters: 1190# gamma - Power (0<gamma<1) of the singularity 1191# f - Function to integrate. Must be a single Tcl command, 1192# to which will be appended the abscissa at which the function 1193# should be evaluated. f is expected to have an integrable 1194# power law singularity at the upper endpoint; that is, the 1195# integrand is expected to diverge as (b-x)**gamma. 1196# a - One bound of the interval 1197# b - The other bound of the interval. a and b need not be in 1198# ascending order. 1199# 1200# Options: 1201# -abserror ABSERROR 1202# Requests that the integration be performed to make 1203# the estimated absolute error of the integral less than 1204# the given value. Default is 1.e-10. 1205# -relerror RELERROR 1206# Requests that the integration be performed to make 1207# the estimated absolute error of the integral less than 1208# the given value. Default is 1.e-6. 1209# -degree N 1210# Specifies the degree of the polynomial that will be 1211# used to extrapolate to a zero step size. -degree 0 1212# requests integration with the midpoint rule; -degree 1 1213# is equivalent to Simpson's 3/8 rule; higher degrees 1214# are difficult to describe but (within reason) give 1215# faster convergence for smooth functions. Default is 1216# -degree 4. 1217# -maxiter N 1218# Specifies the maximum number of triplings of the 1219# number of steps to take in integration. At most 1220# 3**N function evaluations will be performed in 1221# integrating with -maxiter N. The integration 1222# will terminate at that time, even if the result 1223# satisfies neither the -relerror nor -abserror tests. 1224# 1225# Results: 1226# Returns a two-element list. The first element is the estimated 1227# value of the integral; the second is the estimated absolute 1228# error of the value. 1229# 1230# The 'romberg_PowerLawUpper' procedure is a wrapper for 'romberg' 1231# for integrating a function with an integrable power law singularity 1232# at the upper bound of the interval. It works by making the change 1233# of variable u = (b-x)**(1-gamma). 1234# 1235#---------------------------------------------------------------------- 1236 1237proc ::math::calculus::romberg_powerLawUpper { gamma f a b args } { 1238 if { ![string is double -strict $gamma] } { 1239 return -code error [expectDouble $gamma] 1240 } 1241 if { $gamma <= 0.0 || $gamma >= 1.0 } { 1242 return -code error -errorcode [list romberg gammaTooBig] \ 1243 "gamma must lie in the interval (0,1)" 1244 } 1245 if { ![string is double -strict $a] } { 1246 return -code error [expectDouble $a] 1247 } 1248 if { ![string is double -strict $b] } { 1249 return -code error [expectDouble $b] 1250 } 1251 if { $a >= $b } { 1252 return -code error "limits of integration out of order" 1253 } 1254 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 1255 set onemgamma [expr { 1. - $gamma }] 1256 set f [list u_powerLawUpper $f \ 1257 [expr { $gamma / $onemgamma }] \ 1258 [expr { 1. / $onemgamma }] \ 1259 $b] 1260 1261 set limit [expr { pow( $b - $a, $onemgamma ) }] 1262 set result {} 1263 foreach v [eval [linsert $args 0 romberg $f 0 $limit]] { 1264 lappend result [expr { $v / $onemgamma }] 1265 } 1266 return $result 1267 1268} 1269 1270#---------------------------------------------------------------------- 1271# 1272# u_expUpper -- 1273# 1274# Change of variable to integrate a function that decays 1275# exponentially. 1276# 1277# Parameters: 1278# f - Function to integrate 1279# u - Changed variable u = exp(-x) 1280# 1281# Results: 1282# Returns (1/u)*f(-log(u)) 1283# 1284# Side effects: 1285# Whatever f does. 1286# 1287#---------------------------------------------------------------------- 1288 1289proc ::math::calculus::u_expUpper { f u } { 1290 set cmd $f 1291 lappend cmd [expr { -log($u) }] 1292 set y [eval $cmd] 1293 return [expr { $y / $u }] 1294} 1295 1296#---------------------------------------------------------------------- 1297# 1298# romberg_expUpper -- 1299# 1300# Integrate a function that decays exponentially over a 1301# half-infinite interval. 1302# 1303# Parameters: 1304# Same as romberg. The upper limit of integration, 'b', 1305# is expected to be very large. 1306# 1307# Results: 1308# Same as romberg. 1309# 1310# The romberg_expUpper function operates by making the change of 1311# variable, u = exp(-x). 1312# 1313#---------------------------------------------------------------------- 1314 1315proc ::math::calculus::romberg_expUpper { f a b args } { 1316 if { ![string is double -strict $a] } { 1317 return -code error [expectDouble $a] 1318 } 1319 if { ![string is double -strict $b] } { 1320 return -code error [expectDouble $b] 1321 } 1322 if { $a >= $b } { 1323 return -code error "limits of integration out of order" 1324 } 1325 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 1326 set f [list u_expUpper $f] 1327 return [eval [linsert $args 0 \ 1328 romberg $f [expr {exp(-$b)}] [expr {exp(-$a)}]]] 1329} 1330 1331#---------------------------------------------------------------------- 1332# 1333# u_expLower -- 1334# 1335# Change of variable to integrate a function that grows 1336# exponentially. 1337# 1338# Parameters: 1339# f - Function to integrate 1340# u - Changed variable u = exp(x) 1341# 1342# Results: 1343# Returns (1/u)*f(log(u)) 1344# 1345# Side effects: 1346# Whatever f does. 1347# 1348#---------------------------------------------------------------------- 1349 1350proc ::math::calculus::u_expLower { f u } { 1351 set cmd $f 1352 lappend cmd [expr { log($u) }] 1353 set y [eval $cmd] 1354 return [expr { $y / $u }] 1355} 1356 1357#---------------------------------------------------------------------- 1358# 1359# romberg_expLower -- 1360# 1361# Integrate a function that grows exponentially over a 1362# half-infinite interval. 1363# 1364# Parameters: 1365# Same as romberg. The lower limit of integration, 'a', 1366# is expected to be very large and negative. 1367# 1368# Results: 1369# Same as romberg. 1370# 1371# The romberg_expUpper function operates by making the change of 1372# variable, u = exp(x). 1373# 1374#---------------------------------------------------------------------- 1375 1376proc ::math::calculus::romberg_expLower { f a b args } { 1377 if { ![string is double -strict $a] } { 1378 return -code error [expectDouble $a] 1379 } 1380 if { ![string is double -strict $b] } { 1381 return -code error [expectDouble $b] 1382 } 1383 if { $a >= $b } { 1384 return -code error "limits of integration out of order" 1385 } 1386 set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]] 1387 set f [list u_expLower $f] 1388 return [eval [linsert $args 0 \ 1389 romberg $f [expr {exp($a)}] [expr {exp($b)}]]] 1390} 1391 1392 1393# regula_falsi -- 1394# Compute the zero of a function via regula falsi 1395# Arguments: 1396# f Name of the procedure/command that evaluates the function 1397# xb Start of the interval that brackets the zero 1398# xe End of the interval that brackets the zero 1399# eps Relative error that is allowed (default: 1.0e-4) 1400# Result: 1401# Estimate of the zero, such that the estimated (!) 1402# error < eps * abs(xe-xb) 1403# Note: 1404# f(xb)*f(xe) must be negative and eps must be positive 1405# 1406proc ::math::calculus::regula_falsi { f xb xe {eps 1.0e-4} } { 1407 if { $eps <= 0.0 } { 1408 return -code error "Relative error must be positive" 1409 } 1410 1411 set fb [$f $xb] 1412 set fe [$f $xe] 1413 1414 if { $fb * $fe > 0.0 } { 1415 return -code error "Interval must be chosen such that the \ 1416function has a different sign at the beginning than at the end" 1417 } 1418 1419 set max_error [expr {$eps * abs($xe-$xb)}] 1420 set interval [expr {abs($xe-$xb)}] 1421 1422 while { $interval > $max_error } { 1423 set coeff [expr {($fe-$fb)/($xe-$xb)}] 1424 set xi [expr {$xb-$fb/$coeff}] 1425 set fi [$f $xi] 1426 1427 if { $fi == 0.0 } { 1428 break 1429 } 1430 set diff1 [expr {abs($xe-$xi)}] 1431 set diff2 [expr {abs($xb-$xi)}] 1432 if { $diff1 > $diff2 } { 1433 set interval $diff2 1434 } else { 1435 set interval $diff1 1436 } 1437 1438 if { $fb*$fi < 0.0 } { 1439 set xe $xi 1440 set fe $fi 1441 } else { 1442 set xb $xi 1443 set fb $fi 1444 } 1445 } 1446 1447 return $xi 1448} 1449