| … | |
… | |
| 467 | |
467 | |
| 468 | =head3 Inverse Romberg |
468 | =head3 Inverse Romberg |
| 469 | |
469 | |
| 470 | =pod |
470 | =pod |
| 471 | |
471 | |
| 472 | Ueseage: inv_romberg(function_reference, a, value); |
472 | Usage: inv_romberg(function_reference, a, value); |
| 473 | |
473 | |
| 474 | Finds b such that integral from a to b is equal to value. Assumes that the function is positive and continuous. |
474 | Finds b such that the integral of the function from a to b is equal to value. |
|
|
475 | Assumes that the function is continuous and doesn't take on the zero value. |
| 475 | Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals. |
476 | Uses Newton's method of approximating roots of equations, and Romberg to evaluate definite integrals. |
| 476 | |
477 | |
| 477 | =cut |
478 | =cut |
| 478 | |
479 | |
| 479 | sub inv_romberg { |
480 | sub inv_romberg { |
| … | |
… | |
| 484 | my $delta = 1; |
485 | my $delta = 1; |
| 485 | my $i=0; |
486 | my $i=0; |
| 486 | my $funct; |
487 | my $funct; |
| 487 | my $deriv; |
488 | my $deriv; |
| 488 | while (abs($delta) > 0.000001 && $i < 5000) { |
489 | while (abs($delta) > 0.000001 && $i < 5000) { |
| 489 | $funct = romberg($fn_ref,$a,$b)-$value; |
490 | $funct = romberg($fn_ref,$a,$b)-$value; |
| 490 | $deriv = &$fn_ref ( $b ); |
491 | $deriv = &$fn_ref ( $b ); |
|
|
492 | if ($deriv == 0) { |
|
|
493 | warn 'Values of the function are too close to 0.'; |
|
|
494 | return; |
|
|
495 | } |
| 491 | $delta = $funct/$value; |
496 | $delta = $funct/$deriv; |
| 492 | $b = $b - $delta; |
497 | $b = $b - $delta; |
| 493 | $i++; |
498 | $i++; |
| 494 | } |
499 | } |
| 495 | if ($i == 5000) { |
500 | if ($i == 5000) { |
| 496 | warn 'Newtons method does not converge.'; |
501 | warn 'Newtons method does not converge.'; |
| 497 | return; |
502 | return; |
| 498 | } |
503 | } |
| 499 | $b; |
504 | $b; |
| 500 | } |
505 | } |
| 501 | |
506 | |
| 502 | ######################################################### |
507 | ######################################################### |
Parent Directory
| 
Revision Log
| 
Patch