Romberg integration in SciPy - Step-by-Step Execution

Concept Flow - Romberg integration

Define function f(x)

↓

Set integration limits a, b

↓

Initialize Romberg matrix R

↓

Compute R[0,0

↓

For each level k=1 to n

↓

Compute trapezoidal estimates with 2^k intervals

↓

Fill Romberg matrix using Richardson extrapolation

↓

Return R[n,n

↓

END

Romberg integration refines trapezoidal estimates using extrapolation to improve accuracy step-by-step.

Execution Sample

SciPy

import numpy as np
from scipy.integrate import romberg

def f(x):
    return np.sin(x)

result = romberg(f, 0, np.pi, divmax=3)

This code computes the integral of sin(x) from 0 to pi using Romberg integration with up to 3 refinement levels.

Execution Table

Step	k (level)	Number of intervals	Trapezoidal estimate R[k,0]	Richardson extrapolation R[k,j]	Current best estimate
1	0	1	1.5707963267948966	N/A	1.5707963267948966
2	1	2	1.8961188979370398	2.004559754984421	2.004559754984421
3	2	4	1.9742316019455502	2.0002691699483877	2.0002691699483877
4	2	4	N/A	2.0002691699483877	2.0002691699483877
5	3	8	1.999367984224174	2.0000002691699485	2.0000002691699485
6	3	8	N/A	2.0000000000000004	2.0000000000000004
7	3	8	N/A	2.0000000000000004	2.0000000000000004

💡 Reached maximum refinement level divmax=3; final Romberg estimate is R[3,3].

Variable Tracker

Variable	Start	After 1	After 2	After 3	Final
k	0	1	2	3	3
Number of intervals	1	2	4	8	8
R[k,0] (trapezoidal)	1.5707963267948966	1.8961188979370398	1.9742316019455502	1.999367984224174	1.999367984224174
R[k,j] (extrapolated)	N/A	2.004559754984421	2.0002691699483877	2.0000002691699485	2.0000000000000004
Current best estimate	1.5707963267948966	2.004559754984421	2.0002691699483877	2.0000000000000004	2.0000000000000004

Key Moments - 3 Insights

Why does the Romberg matrix have entries labeled R[k,j] and what do they represent?

R[k,0] is the trapezoidal estimate with 2^k intervals; R[k,j] for j>0 are refined estimates using Richardson extrapolation combining previous values, as shown in execution_table rows 2-7.

Why does the number of intervals double at each level k?

Each level k refines the integral by doubling intervals to get a better trapezoidal estimate R[k,0], visible in the 'Number of intervals' column increasing 1,2,4,8.

Why does the final estimate improve and approach the true integral value?

Because Richardson extrapolation reduces error by combining trapezoidal estimates, the final R[k,k] converges to the true integral, as seen by the 'Current best estimate' approaching 2.

Visual Quiz - 3 Questions

Look at the execution_table at step 3. What is the trapezoidal estimate R[2,0]?

A 1.9742316019455502

B 1.8961188979370398

C 2.004559754984421

D 1.5707963267948966

💡 Hint

Check the 'Trapezoidal estimate R[k,0]' column at step 3 where k=2.

At which step does the Romberg estimate first use 8 intervals?

A Step 4

B Step 5

C Step 3

D Step 2

💡 Hint

Look at the 'Number of intervals' column to find when it reaches 8.

If divmax was set to 2 instead of 3, what would be the final Romberg estimate step?

A Step 3

B Step 5

C Step 4

D Step 2

💡 Hint

The final step corresponds to k=divmax; check the 'k (level)' column.

Concept Snapshot

Romberg integration:
- Uses trapezoidal rule estimates with 2^k intervals
- Builds a Romberg matrix R[k,j] using Richardson extrapolation
- R[k,0] is trapezoidal estimate; R[k,j] refines it
- Final estimate is R[n,n] for n=divmax
- Improves accuracy by canceling error terms
- scipy.integrate.romberg automates this process

Full Transcript

Romberg integration is a method to calculate integrals more accurately by combining trapezoidal rule estimates with increasing numbers of intervals. We start by defining the function to integrate and the limits. Then, we compute trapezoidal estimates with 1, 2, 4, 8 intervals, etc. These estimates fill the first column of a Romberg matrix. Using Richardson extrapolation, we fill the rest of the matrix to refine the integral estimate. Each step improves accuracy by reducing error. The final value at the bottom-right of the matrix is the best estimate. This process is automated in scipy's romberg function. The execution table shows each step's trapezoidal estimate, extrapolated values, and the current best estimate. Variables like k (level) and number of intervals double each step. Beginners often wonder about the meaning of R[k,j], why intervals double, and how accuracy improves. The quiz questions help check understanding of these steps.