Exercises — Implementing root-finding from scratch — Newton-Raphson, bisection

Humein kya chahiye: ek bracket $[a,b]$ tabhi valid hai jab $==f(a)\cdot f(b)<0==$ — function sign change karta hai, toh Intermediate Value Theorem ke hisaab se ek continuous curve ko andar zero cross karna hi padega.

Endpoints evaluate karo:

• $f(0)=-2$, $f(1)=1-1-2=-2$ → product $=+4>0$. Valid nahi (dono negative — curve yahan zero tak dip nahi bhi kar sakti).
• $f(1)=-2$, $f(2)=8-2-2=4$ → product $=-8<0$. Valid ✓ (negative → positive, ek guaranteed crossing).
• $f(2)=4$, $f(3)=27-3-2=22$ → product $=+88>0$. Valid nahi.

Answer: sirf $[1,2]$ ek root bracket karta hai.

Neeche di gayi figure curve ko is valid bracket ke saath dikhati hai, iske pehle bisection midpoint (plum color mein), aur Exercise 2.2 ka Newton tangent (orange) — abhi ek baar dekho ki dono methods ek hi crossing ko kaise alag tarike se approach karte hain.

Dekho ki plum midpoint $m=1.5$ par pehle se zero ke kaafi kareeb hai (kyunki curve yahan almost seedhi hai), jabki orange tangent ka x-intercept $x_0=2$ se ek hi jump mein true root ke almost upar aa jaata hai.

• (a) Bisection. Newton ko derivative chahiye; bisection ko sirf $f$ aur ek sign change chahiye.
• (b) Newton. Close start + derivative = quadratic convergence, kam steps.
• (c) Bisection. Tiny $f^{\prime }(x_n)$ Newton step $-f(x_n)/f^{\prime }(x_n)$ ko explode kara deta hai — tangent almost flat hai aur iska x-intercept bahut door ud jaata hai.

Store karo $f_a=f(1)=-2$ (negative), $f_b=f(2)=4$ (positive).

Step 1: midpoint $m=(1+2)/2=1.5$. Compute karo $f(1.5)=3.375-1.5-2=-0.125$ (negative). Stored $f_a$ se compare karo: $f_a\cdot f_m=(-2)(-0.125)=+0.25>0$ → left side par same sign, toh crossing right half mein hai. Naya bracket: $[1.5,2]$, aur ab $f_a=f(1.5)=-0.125$.

Step 2: midpoint $m=(1.5+2)/2=1.75$. Compute karo $f(1.75)=5.359375-1.75-2=1.609375$ (positive). $f_a\cdot f_m=(-0.125)(1.609)<0$ → left par sign change, left half rakho. Naya bracket: $[1.5,1.75]$.

Agar tum abhi rokte to midpoint $(1.5+1.75)/2=1.625$ return karte. (True root $\approx 1.5214$ hai.)

Update rule: $x_{n+1}=x_n-\frac{x_n^3-x_n-2}{3x_n^2-1}$.

Step 1 ($x_0=2$): $f(2)=4$, $f^{\prime }(2)=3\cdot 4-1=11$. $x_1=2-\frac{4}{11}=2-0.363636\dots =\mathrm{1.636364.}$

Step 2 ($x_1=1.636364$): $f(x_1)=1.636364^3-1.636364-2=4.380697-3.636364=0.744333$; $f^{\prime }(x_1)=3(1.636364{)}^2-1=8.032231-1=7.032231$. $x_2=1.636364-\frac{0.744333}{7.032231}=1.636364-0.105844=\mathrm{1.530520.}$

Do Newton steps ne humein $1.53052$ tak pahuncha diya (true root $1.52138$); do bisection steps sirf $1.625$ tak pahunche. Newton fast aage nikal jaata hai jab guess kareeb ho.

Bound derive karna (geometric "kyun"): Interval width $w_0=b-a$ se shuru karo. Har bisection step live interval ko exactly aadha kar deta hai, kyunki hum midpoint par split karte hain aur ek side rakhte hain. Toh $n$ steps ke baad width hai: $w_n=\frac{w_0}{2^n}=\frac{b-a}{2^n}.$ Ab — us final interval ke andar true root kahan hai? Hum surviving interval ka midpoint return karte hain. Root kahi bhi interval mein ho sakta hai, aur midpoint dead-centre baitha hai, toh root midpoint se sabse zyada interval width ka aadha ho sakta hai (root ek end par, midpoint beech mein). Isliye error mein extra factor of two hai: $\text{error}\le \frac{w_n}{2}=\frac{b-a}{2^{n+1}}.$ $2^{n+1}$ koi magic constant nahi hai: $2^n$ halvings count karta hai, aur extra $2$ hai "midpoint half-width se zyada door nahi" geometry.

Step count ke liye solve karna. Humein chahiye: $\frac{b-a}{2^{n+1}}\le \epsilon ⟺\frac{1}{2^{n+1}}\le 10^{-6}⟺2^{n+1}\ge 10^6.$ ${\log }_2$ lo: $n+1\ge {\log }_2(10^6)=6{\log }_{2}10\approx 6\times 3.32193=19.9316$. Toh $n\ge 18.93$, matlab $n=19$ steps kaafi hain. (Har step ek bit buy karta hai; $2^{20}=1048576>10^6$.)

Taylor's theorem se $C$ derive karna. Maano $e_n=x_n-x^{\ast }$ error hai, jahan $x^{\ast }$ true root hai ($f(x^{\ast })=0$). $f$ aur $f^{\prime }$ ko $x^{\ast }$ ke aaspaas Taylor-expand karo, $e_n^2$ tak terms rakhte hue: $f(x_n)=f(x^{\ast })+f^{\prime }(x^{\ast })e_n+\frac{1}{2}{f}^{\prime \prime }(x^{\ast })e_n^2+\cdots =f^{\prime }(x^{\ast })e_n+\frac{1}{2}{f}^{\prime \prime }(x^{\ast })e_n^2,$ ($f(x^{\ast })$ term vanish ho jaata hai kyunki $x^{\ast }$ root hai), aur $f^{\prime }(x_n)=f^{\prime }(x^{\ast })+f^{\prime \prime }(x^{\ast })e_n+\cdots .$ Agla Newton error hai $e_{n+1}=x_{n+1}-x^{\ast }=e_n-\frac{f(x_n)}{f^{\prime }(x_n)}$. Dono expansions substitute karo: $e_{n+1}=e_n-\frac{f^{\prime }(x^{\ast })e_n+\frac{1}{2}{f}^{\prime \prime }(x^{\ast })e_n^2}{f^{\prime }(x^{\ast })+f^{\prime \prime }(x^{\ast })e_n}.$ Numerator se $f^{\prime }(x^{\ast })e_n$ aur denominator se $f^{\prime }(x^{\ast })$ factor out karo: $e_{n+1}=e_n-e_n\cdot \frac{1+\frac{f^{\prime \prime }(x^{\ast })}{2f^{\prime }(x^{\ast })}{e}_n}{1+\frac{f^{\prime \prime }(x^{\ast })}{f^{\prime }(x^{\ast })}{e}_n}.$ Ab chhote $e_n$ ke liye $\frac{1+\alpha e_n}{1+\beta e_n}\approx (1+\alpha e_n)(1-\beta e_n)\approx 1+(\alpha -\beta )e_n$ use karo (denominator ka geometric-series approximation). $\alpha =\frac{f^{\prime \prime }}{2f^{\prime }}$ aur $\beta =\frac{f^{\prime \prime }}{f^{\prime }}$ ke saath: $e_{n+1}\approx e_n-e_n(1+(\alpha -\beta )e_n)=-e_n^2(\alpha -\beta )=e_n^2(\beta -\alpha ).$ Aur $\beta -\alpha =\frac{f^{\prime \prime }}{f^{\prime }}-\frac{f^{\prime \prime }}{2f^{\prime }}=\frac{f^{\prime \prime }}{2f^{\prime }}$. Toh constant hai: $\boxed{e_{n+1}\approx \frac{f^{\prime \prime }(x^{\ast })}{2f^{\prime }(x^{\ast })}{e}_n^2}\Rightarrow C=\frac{f^{\prime \prime }(x^{\ast })}{2f^{\prime }(x^{\ast })}.$ Isliye $C$ ka exactly wahi form hai — yeh leading curvature term $f^{\prime \prime }/2$ aur slope $f^{\prime }$ ka ratio hai, seedha Taylor remainder se nikla hua.

Yahan apply karna. $f(x)=x^2-2$ ke liye: $f^{\prime }(x)=2x$ toh $f^{\prime }(x^{\ast })=2\sqrt{2}$; $f^{\prime \prime }(x)=2$ toh $f^{\prime \prime }(x^{\ast })=2$. $C=\frac{2}{2\cdot 2\sqrt{2}}=\frac{1}{2\sqrt{2}}\approx \mathrm{0.353553.}$ Predicted $e_1\approx C{e}_0^2=0.353553\times (0.085786{)}^2=0.353553\times 0.0073593=0.0026019$.

Actual: $x_1=\frac{1.5}{2}+\frac{1}{1.5}=0.75+0.666667=1.416667$, toh $e_1=1.416667-1.414214=0.002453$.

Predicted $0.002602$ vs actual $0.002453$ — same order, ~$6\%$ apart (difference hai woh higher-order terms jo humne expansion mein drop kiye). Prediction squaring behaviour confirm karta hai: $e_0\approx 0.086$ se $e_1\approx 0.0025$ ho gaya, roughly $e_0^2$ scaled — accuracy ke digits double ho gaye.

def safe_step(f, df, a, b, x, fa):
    dfx = df(x)
    if dfx != 0:
        xn = x - f(x) / dfx        # try Newton
        if a < xn < b:             # accept only if it stays inside the bracket
            x_new = xn
        else:
            x_new = (a + b) / 2    # Newton escaped -> bisect
    else:
        x_new = (a + b) / 2        # flat tangent -> bisect
    # shrink the bracket using the new point's sign
    fx = f(x_new)
    if fa * fx < 0:
        b = x_new                  # sign change on the left half
    else:
        a, fa = x_new, fx          # sign change on the right half
    return a, b, x_new, fa

Har guard kyun hai:

• dfx != 0: flat tangent infinite/undefined step deta hai — Newton ka exact failure mode. Bisection par fall back karo, jo kabhi slope se divide nahi karta.
• a < xn < b: agar Newton ka tangent trapped region ke bahar point kare, toh root wahan ho nahi sakta — guess ek blind midpoint se bhi worse hai. Reject karo.
• Sign test se a/b update karna valid bracket hamesha maintain karta hai, toh hum bisection ki guarantee inherit karte hain aur Newton ki speed bhi enjoy karte hain jab woh behave kare.

Hand trace $f(x)=x^3-x-2$ par, $[1,2]$, start $x=2$, $f_a=f(1)=-2$:

• $df(2)=11$, Newton $x_n=2-4/11=1.6364$, jo $[1,2]$ ke andar hai → accept. $f(1.6364)=0.7443>0$; $f_a\cdot f_x=(-2)(0.744)<0$ → right end shrink karo: bracket $[1,1.6364]$ ho jaata hai.

Yahi philosophy hai scipy.optimize.brentq ke peeche: bracket rakho, andar fast interpolation use karo, jab fast step misbehave kare toh bisect karo.

Newton hai $x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}$. Jab $f^{\prime }$ ka formula na ho, slope ko last do points ke through line se approximate karo (ek finite difference, derivative ki definition se closely related): $f^{\prime }(x_n)\approx \frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}.$ Substitute karo: $x_{n+1}=x_n-f(x_n)\cdot \frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}.$ Gain: koi derivative nahi chahiye, aur sirf ek new $f$-evaluation per step (Newton ko $f$ aur $f^{\prime }$ dono chahiye). Cost: convergence order $2$ se golden ratio $\phi \approx 1.618$ tak drop ho jaata hai — phir bhi superlinear, bisection se faster, true Newton se slower. Yeh Newton ki fragility bhi inherit karta hai (koi bracket guarantee nahi).

Edge case — vanishing denominator. Formula $f(x_n)-f(x_{n-1})$ se divide karta hai. Agar do most recent function values equal (ya nearly so) hon, toh yeh denominator zero (ya tiny) hai aur step blow up karta hai ya division error ke saath crash karta hai. Geometrically yeh ek horizontal secant line hai: ek flat line kabhi x-axis cross nahi karti, toh "secant zero kahan hit karta hai?" ka koi finite answer nahi — Newton ke flat-tangent failure ka exact analogue. Yeh ek local extremum ke kareeb hota hai, ya jab do iterates accidentally same height par land kar jaate hain.

Guard:

def secant_step(f, x_prev, x_cur, f_prev, f_cur, a, b):
    denom = f_cur - f_prev
    if denom == 0:                        # horizontal secant -> no intercept
        return (a + b) / 2                # fall back to a bracketed bisection
    return x_cur - f_cur * (x_cur - x_prev) / denom

Mastery level par Exercise 4.1 ka wahi lesson apply hota hai: raw secant fragile hai, toh production code ise safeguard karta hai — denominator check karo, aur bracket rakho jab secant point escape kare ya denominator collapse kare. Yahi denominator guard plus bracket fallback exactly wajah hai ki real solvers jaise scipy.optimize.brentq kabhi degenerate slope se divide nahi karte.

Yahan $f(x)=x^{1/3}$ hai toh $f^{\prime }(x)=\frac{1}{3}{x}^{-2/3}$. Newton step: $x_{n+1}=x_n-\frac{x_n^{1/3}}{\frac{1}{3}{x}_n^{-2/3}}=x_n-3x_n^{1/3}{x}_n^{2/3}=x_n-3x_n=-2x_n.$ Toh har iteration $x_{n+1}=-2x_n$ deta hai: guess magnitude mein double ho jaata hai aur sign flip karta hai har step par. $x_0=1$ se: $1\to -2\to 4\to -8\to \dots$ — woh $0$ se hamesha door bhaag jaata hai.

Geometry: $f$ ka root par infinitely steep tangent hai (vertical inflection). $0$ se door tangent itna shallow hai ki iska x-intercept jahan tum shuru hue wahan se zyada door land karta hai. Neeche di gayi figure yeh runaway sketch karti hai — teen dashed tangents (orange, plum, teal) follow karo: har ek ka x-intercept (triangle marker) us point se zyada door baitha hai jahan se woh aaya.

Figure ko left-to-right padho: $x_0=1$ se $-2$ par intercept milta hai, jiske tangent se $4$ par intercept milta hai, aur aise hi — markers dono directions mein baahir jaate hain, kabhi central star $0$ ki taraf nahi.

Fix: Bisection use karo. Kisi bhi bracket jaise $[-1,2]$ (signs $-$ aur $+$) ke saath, yeh inexorably $0$ ki taraf march karta hai. Yahi parent note ki headline warning hai: bisection guaranteed hai, Newton nahi.

Parent code yaad karo:

if fa == 0:
    return a
...
if fa * fb > 0:
    raise ValueError(...)

$f(a)=0$ ke saath: pehla guard fa == 0 fire karta hai aur a immediately return karta hai — correct, yeh exact root hai.

Order kyun matter karta hai: agar humne pehle fa * fb > 0 check kiya hota, note karo $f(a)\cdot f(b)=0\cdot f(b)=0$, jo na $<0$ hai na $>0$. Yeh strict bracket test fail karta hai (> 0 False hai, toh koi error raise nahi hota) lekin loop ka termination (b-a)/2 < tol true root ka return kabhi trigger nahi kar sakta, aur isse bhi bura, sign-update logic genuine strict sign change assume karta hai. Ek zero endpoint us assumption ko break karta hai. Endpoints pehle check karna poore degenerate case ko cleanly sidestep karta hai. Yeh exactly parent note mein second steel-man hai.

Step $10^{-11}$ bada hai $\text{tol}=10^{-12}$ se, toh abs(step) < tol test False hai → loop ek aur iteration continue karta hai. (Agla step ke baad yeh almost certainly tol se neeche aa jaayega aur ruk jaayega.)

$\text{tol}=10^{-20}$ push karne par: reliably meet karna impossible. Double-precision floats lagbhag $15$–$16$ significant decimal digits carry karte hain; $x\approx 1.41$ ke kareeb representable numbers ke beech spacing (machine epsilon scale, value ka $\approx 2.2\times 10^{-16}$ times) roughly $2\times 10^{-16}$ hai. Tum us gap se chhota step resolve nahi kar sakte — loop maxit tak spin karta rahega bina kabhi $|step|<10^{-20}$ satisfy kiye. Tolerances machine epsilon ke kuch times ya usse upar set karo; number system se zyada precision maangna ek guaranteed non-termination bug hai.