|
Last change
on this file since 7837 was
7837,
checked in by mungkasi, 14 years ago
|
|
Again, adding some codes for 1d problems on debris avalanche and periodic waves.
|
|
File size:
902 bytes
|
| Line | |
|---|
| 1 | ## module bisect |
|---|
| 2 | ''' root = bisect(f,x1,x2,switch=0,tol=1.0e-9). |
|---|
| 3 | Finds a root of f(x) = 0 by bisection. |
|---|
| 4 | The root must be bracketed in (x1,x2). |
|---|
| 5 | Setting switch = 1 returns root = None if |
|---|
| 6 | f(x) increases as a result of a bisection. |
|---|
| 7 | ''' |
|---|
| 8 | from math import log,ceil |
|---|
| 9 | import error |
|---|
| 10 | |
|---|
| 11 | def bisect(f,x1,x2,switch=0,epsilon=1.0e-15):#9): |
|---|
| 12 | f1 = f(x1) |
|---|
| 13 | if f1 == 0.0: return x1 |
|---|
| 14 | f2 = f(x2) |
|---|
| 15 | if f2 == 0.0: return x2 |
|---|
| 16 | if f1*f2 > 0.0: error.err('Root is not bracketed') |
|---|
| 17 | n = ceil(log(abs(x2 - x1)/epsilon)/log(2.0)) |
|---|
| 18 | for i in range(n): |
|---|
| 19 | x3 = 0.5*(x1 + x2); f3 = f(x3) |
|---|
| 20 | if (switch == 1) and (abs(f3) >abs(f1)) \ |
|---|
| 21 | and (abs(f3) > abs(f2)): |
|---|
| 22 | return None |
|---|
| 23 | if f3 == 0.0: return x3 |
|---|
| 24 | if f2*f3 < 0.0: |
|---|
| 25 | x1 = x3; f1 = f3 |
|---|
| 26 | else: |
|---|
| 27 | x2 = x3; f2 = f3 |
|---|
| 28 | return (x1 + x2)/2.0 |
|---|
Note: See
TracBrowser
for help on using the repository browser.
