scipy.integrate — odeint, solve_ivp (RK45, DOP853), quad

5.4.9 · Coding › Scientific Computing (Python)

1. `quad` — numerical definite integral

Simplest rule ko scratch se derive karna (trapezoid → Simpson)

Hum chahte hain $\int_{a}^{b} f d x$ sirf $f (a), f (m), f (b)$ use karke, jahan $m = \frac{a + b}{2}$ . In 3 points se ek parabola $p (x) = α + β x + γ x^{2}$ fit karo aur use exactly integrate karo. Shift karo taaki $a = - h, m = 0, b = h$ . Phir:

$\int_{- h}^{h} (α + β x + γ x^{2}) d x = 2 α h + \frac{2}{3} γ h^{3}$

Odd term vanish ho jaata hai (Why? $\int_{- h}^{h} x d x = 0$ symmetry se). Ab $α, γ$ ko samples ke through express karo: $f (0) = α$ , aur $f (h) + f (- h) = 2 α + 2 γ h^{2} \Rightarrow γ h^{2} = \frac{f ( h ) + f ( - h )}{2} - f (0)$ . Substitute karo:

$\int_{- h}^{h} f \approx \frac{h}{3} [f (- h) + 4 f (0) + f (h)]$

Yahi hai Simpson's rule. quad yahi idea hai, lekin adaptive aur higher order.

from scipy.integrate import quad
import numpy as np
 
val, err = quad(lambda x: np.exp(-x**2), 0, np.inf)
# val ≈ 0.8862269  (= sqrt(pi)/2), err ≈ 7.1e-09

val = estimate, err = estimated absolute error.
Extra params args=(...) se pass karo; kinks flag karne ke liye points=[...] use karo.

2. `solve_ivp` — modern ODE solver

Ek single step derive karna (forward Euler, phir RK)

Exact statement (calculus ka fundamental theorem) se shuru karo: $y (t_{n} + h) = y (t_{n}) + \int_{t_{n}}^{t_{n} + h} f (t, y) d t$

WHY yeh seed hai: agar hum us integral ko exactly jaante toh kaam ho jaata; har method bas us chhote integral ke liye ek quadrature rule hai.

Euler (left-rectangle): integral ko $h f (t_{n}, y_{n})$ se approximate karo: $y_{n + 1} = y_{n} + h f (t_{n}, y_{n}) (error O (h^{2}) per step)$
RK4 / RK45 (Runge–Kutta): step ke andar kuch trial points pe slope sample karo, jaise Simpson ne ek interval ke andar kiya, aur blend karo: $y_{n + 1} = y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$ jahan $k_{1} = f (t_{n}, y_{n})$ , $k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1})$ , wagera.

Worked example Stiff problems aur method menu

method	order	adaptive	use when
`RK45`	5(4)	yes	default, smooth non-stiff
`DOP853`	8	yes	high accuracy chahiye
`Radau`, `BDF`, `LSODA`	implicit	yes	stiff systems

Ek system stiff hota hai jab usmein bahut fast aur bahut slow time-scales saath hon; explicit RK stability ke liye absurdly tiny $h$ force karta hai, isliye implicit BDF/Radau use karo.

Worked example — exponential decay $\overset{y}{˙} = - 2 y, y (0) = 3$

Exact solution: $y (t) = 3 e^{- 2 t}$ .

from scipy.integrate import solve_ivp
import numpy as np
 
sol = solve_ivp(lambda t, y: -2*y, [0, 5], [3.0],
                method='RK45', t_eval=np.linspace(0, 5, 100),
                rtol=1e-8, atol=1e-10)
sol.t      # times
sol.y[0]   # y values, shape (n_states, n_times)

Why y0=[3.0] (ek list)? solve_ivp hamesha state ko ek vector maanta hai; scalar bhi 1-element hona chahiye.
Why sol.y[0]? rows state components hain, columns time points hain.

Worked example — ek system (SHM / 2nd-order → 1st-order)

$\overset{x}{¨} = - ω^{2} x$ . Why convert karte hain? Solvers sirf first-order systems lete hain. $y_{0} = x, y_{1} = \overset{x}{˙}$ rakh lo: $\overset{y}{˙}_{0} = y_{1}, \overset{y}{˙}_{1} = - ω^{2} y_{0}$

w = 2.0
def rhs(t, y):
    return [y[1], -w**2 * y[0]]
sol = solve_ivp(rhs, [0, 10], [1, 0], method='DOP853', rtol=1e-10)
# x(t) = sol.y[0] should match cos(2t); energy conserved

3. `odeint` — legacy classic

from scipy.integrate import odeint
y = odeint(lambda y, t: -2*y, 3.0, np.linspace(0,5,100))  # note (y, t)!

Shape (len(t), n_states) ka array return karta hai — yeh bhi solve_ivp se transposed hai.
Modern advice: solve_ivp prefer karo (events, dense output, method choice). odeint ki knowledge purana code padhne ke liye rakhho.

Common mistake Classic errors ko steel-man karna

1. Signature swap karna. Lagta hai sahi hai ki solve_ivp ke liye f(y, t) likho kyunki odeint aisa karta hai, aur "y hi main cheez hai." Fix: solve_ivp mein f(t, y) hai (time pehle); odeint mein f(y, t) hai. Neeche mnemonic hai.

2. "Main default tolerance pe trust kar lunga." Safe lagta hai kyunki call succeed hoti hai. Lekin default rtol=1e-3 oscillators/long times ke liye visibly wrong ho sakta hai. Fix: jab accuracy matter kare toh rtol/atol explicitly set karo.

3. t_eval ko integration grid samajh lena. Lagta hai t_eval step size control karta hai. Nahi karta — solver khud adaptive steps choose karta hai; t_eval sirf choose karta hai ki result kahan report ho. Fix: accuracy badlani ho toh tolerances ya method badlo, t_eval nahi.

4. Stiff problem pe RK45 use karna aur use crawl/blow up hote dekhna. Bug lagta hai. Yeh stability hai: explicit methods ko tiny $h$ chahiye. Fix: BDF/Radau/LSODA pe switch karo.

Recall Feynman: 12-saal ke bacche ko samjhao

Socho tumhe pata hai ki ek toy car har pal kitni tez chal rahi hai, lekin yeh nahi ki woh kahan hai. Yeh pata karne ke liye ki woh kahan pahunchi, ek tiny moment lo, uski current speed se thoda aage badhao, speed phir check karo, aur repeat karo — yahi hai solve_ivp jo ODE solve karta hai. Smart solvers har tiny step ke andar kuch baar peek karte hain (RK45) taaki zyada accurately aage badh sakein, aur bade steps lete hain smooth seedhi sadak pe aur chhote steps curvy jagah pe. quad alag hai: wahan tumhe pehle se ek hill ki shape pata hai aur bas uske neeche ka area chahiye, toh clever jagahon pe height sample karte ho aur add karte ho.

Flashcards

quad kya return karta hai?

Ek tuple (value, estimated_absolute_error).

solve_ivp ke RHS ka argument order?

f(t, y) — time pehle.

odeint ke RHS ka argument order?

f(y, t) — state pehle (solve_ivp se ulta).

RK45 ka default kya hai aur yeh step size kaise choose karta hai?

Dormand–Prince 5(4): ek embedded 4th-order estimate ko 5th-order se compare karta hai; difference local error hai jo

h

adapt karne ke liye use hota hai. :::

Stiff problem pe RK45/DOP853 se BDF/Radau pe kab switch karna chahiye?

Jab ODE stiff ho (bahut alag time-scales) aur explicit methods ko impractically chhote steps chahiye hon. :::

t_eval kya control karta hai (aur kya NAHI karta)?

Kahan solution report ho; yeh solver ke internal adaptive step size ko set NAHI karta. :::

[- h, h]

pe Simpson's rule formula?

\frac{h}{3} [f (- h) + 4 f (0) + f (h)]

, cubics ke liye exact. :::

2nd-order ODE ko solve karne se pehle kyun rewrite karna padta hai?

Solvers sirf first-order systems accept karte hain;

y_{1} = \overset{x}{˙}

introduce karo taaki

\overset{y}{˙}_{0} = y_{1}, \overset{y}{˙}_{1} = \overset{x}{¨}

ho. :::

solve_ivp solution array sol.y ki shape?

(n_states, n_time_points) — components rows hain. :::

DOP853 apni cost ke liye kyun worth hai?

8th order: smooth problems pe bade accurate steps → tight tolerances ke liye kam evaluations. :::

odeint aur solve_ivp dono ka output shape kya hai?

odeint: (len(t), n_states); solve_ivp: (n_states, n_time_points) — dono transposed hain ek doosre se. :::

Simpson's rule mein odd term kyun vanish hota hai?

Kyunki

\int_{- h}^{h} x d x = 0

symmetry se — odd functions ka symmetric interval pe integral zero hota hai. :::

Gauss quadrature

n

points ke saath kitne degree tak exact hai?

Degree

2 n - 1

tak ke polynomials ke liye exact. :::

solve_ivp mein y0 ko list kyun banana padta hai scalar ke liye bhi?

solve_ivp state ko hamesha vector maanta hai; isliye scalar bhi 1-element list/array hona chahiye. :::

Simpson's rule ka error term kis derivative pe depend karta hai aur kyun?

f^{(4)}

pe — kyunki humne degree-2 polynomial match ki, toh pehla uncaptured term 4th derivative (next even derivative) se aata hai. :::

Stiff system ka matlab kya hai?

Woh system jisme bahut fast aur bahut slow time-scales saath hoon, jisse explicit methods stability ke liye absurdly tiny

h

use karne ko majboor hoon. :::

solve_ivp mein t_eval set karne se accuracy kyun nahi badh'ti?

Kyunki internal adaptive steps solver khud choose karta hai tolerances ke basis pe; t_eval sirf reporting points set karta hai. :::

Connections

Riemann sums aur Fundamental Theorem of Calculus — yahan sab integration ka seed.
Runge-Kutta methods · Stiff differential equations · Numerical stability
Gaussian quadrature — quad ke peeche ka math.
numpy.trapz and simpson — jab sirf samples ho toh fixed-grid integration.
Harmonic oscillator — canonical solve_ivp test system.