Simpson revival
08-04-2018, 05:29 PM
Post: #15
 Albert Chan Senior Member Posts: 1,676 Joined: Jul 2018
RE: Simpson revival
In Python, Romberg's method corrections can be simplified with generator.

Code:
from __future__ import division def simpson(f, a, b, n=2):          # Simpson generator     h = (b-a) / 2     t0 = (f(a)+f(b)) * h     while 1:         s = sum(f(a + i*h) for i in xrange(1, n, 2))         t1 = t0/2 + s*h             # simpson's partial sum         yield t1 + (t1-t0)/3        # simpson's formula         t0, n, h = t1, 2*n, h/2

Code:
def _romberg(gen, y, k):            # Romberg generator     for z in gen:         yield z + (z-y)/k         y = z

Code:
def romberg(f, a, b, verbal=False, k=15):     "Romberg's Method integration"     gen = simpson(f, a, b, n=2)     x = gen.next()     while 1:         y = gen.next()         x = y + (y-x) / k         if verbal: print x         if x == y: return x         # stop if converged         gen = _romberg(gen, y, k)   # Romberg corrections         k = 4 * k + 3

Redo previous post example:

>>> romberg(lambda x: 1/x, 1, 2, verbal=True)
0.693174603175
0.693147477645
0.693147181917
0.693147180562
0.69314718056
0.69314718056
0.69314718055994506
 « Next Oldest | Next Newest »

 Messages In This Thread Simpson revival - Pekis - 09-28-2016, 08:04 AM RE: Simpson revival - Namir - 09-28-2016, 10:40 AM RE: Simpson revival - Dieter - 09-28-2016, 12:25 PM RE: Simpson revival - Namir - 09-28-2016, 02:00 PM RE: Simpson revival - Namir - 09-28-2016, 03:50 PM RE: Simpson revival - Dieter - 09-28-2016, 10:00 PM RE: Simpson revival - Albert Chan - 07-31-2018, 02:57 PM RE: Simpson revival - Namir - 09-29-2016, 04:44 PM RE: Simpson revival - Dieter - 09-29-2016, 06:23 PM RE: Simpson revival - Namir - 09-29-2016, 10:27 PM RE: Simpson revival - Dieter - 09-30-2016, 06:00 AM RE: Simpson revival - Dieter - 10-02-2016, 03:29 PM RE: Simpson revival - Namir - 10-02-2016, 04:48 PM RE: Simpson revival - Namir - 10-09-2016, 03:14 AM RE: Simpson revival - Albert Chan - 08-04-2018 05:29 PM

User(s) browsing this thread: 1 Guest(s)