Intuition The big pictureSome ODEs contain very fast decaying parts (e − 1000 t e^{-1000t} e − 1000 t slow parts (e − t e^{-t} e − t explicit method (like forward Euler) is forced to
take microscopic steps — not for accuracy, but just to stay stable . That is what we call a
stiff equation. The cure is implicit methods (backward Euler), which stay stable with big
steps. The whole topic is a battle of stability vs step size .
An ODE is stiff when it has solution components decaying on wildly different time scales ,
so that an explicit method needs a step size h h h fastest (smallest) time scale,
even where the true solution is smooth and slow .
Take the model test problem (the lens through which we study EVERY method):
y ′ = λ y , y ( 0 ) = y 0 , λ < 0 ( or Re λ < 0 ) y' = \lambda y, \qquad y(0)=y_0, \qquad \lambda < 0 \ (\text{or } \operatorname{Re}\lambda<0) y ′ = λ y , y ( 0 ) = y 0 , λ < 0 ( or Re λ < 0 )
True solution: y ( t ) = y 0 e λ t → 0 y(t)=y_0 e^{\lambda t}\to 0 y ( t ) = y 0 e λ t → 0 λ \lambda λ faster the decay,
the stiffer the problem.
Derive forward Euler from the definition of the derivative:
y ′ ( t ) ≈ y n + 1 − y n h ⇒ y n + 1 − y n h = f ( t n , y n ) y'(t)\approx \frac{y_{n+1}-y_n}{h} \;\Rightarrow\; \frac{y_{n+1}-y_n}{h}=f(t_n,y_n) y ′ ( t ) ≈ h y n + 1 − y n ⇒ h y n + 1 − y n = f ( t n , y n )
y n + 1 = y n + h f ( t n , y n ) \boxed{y_{n+1}=y_n+h\,f(t_n,y_n)} y n + 1 = y n + h f ( t n , y n )
Apply to y ′ = λ y y'=\lambda y y ′ = λ y
y n + 1 = y n + h λ y n = ( 1 + h λ ) y n ⇒ y n = ( 1 + h λ ) n y 0 y_{n+1}=y_n+h\lambda y_n=(1+h\lambda)\,y_n \;\Rightarrow\; y_n=(1+h\lambda)^n y_0 y n + 1 = y n + hλ y n = ( 1 + hλ ) y n ⇒ y n = ( 1 + hλ ) n y 0
If λ = − 1000 \lambda=-1000 λ = − 1000 h < 0.002 h<0.002 h < 0.002 forever , even after the transient is gone. Disaster.
Evaluate f f f new point t n + 1 t_{n+1} t n + 1 y n + 1 y_{n+1} y n + 1 both sides — you must solve an equation. That extra work buys huge stability.
Derive backward Euler using a backward difference at t n + 1 t_{n+1} t n + 1
y ′ ( t n + 1 ) ≈ y n + 1 − y n h = f ( t n + 1 , y n + 1 ) y'(t_{n+1})\approx \frac{y_{n+1}-y_n}{h} = f(t_{n+1},y_{n+1}) y ′ ( t n + 1 ) ≈ h y n + 1 − y n = f ( t n + 1 , y n + 1 )
y n + 1 = y n + h f ( t n + 1 , y n + 1 ) \boxed{y_{n+1}=y_n+h\,f(t_{n+1},\,y_{n+1})} y n + 1 = y n + h f ( t n + 1 , y n + 1 )
This is implicit : y n + 1 y_{n+1} y n + 1 y ′ = λ y y'=\lambda y y ′ = λ y
y n + 1 = y n + h λ y n + 1 ⇒ ( 1 − h λ ) y n + 1 = y n y_{n+1}=y_n+h\lambda y_{n+1}\;\Rightarrow\;(1-h\lambda)y_{n+1}=y_n y n + 1 = y n + hλ y n + 1 ⇒ ( 1 − hλ ) y n + 1 = y n
y n + 1 = 1 1 − h λ y n ⇒ y n = ( 1 1 − h λ ) n y 0 y_{n+1}=\frac{1}{1-h\lambda}\,y_n \;\Rightarrow\; y_n=\Big(\frac{1}{1-h\lambda}\Big)^n y_0 y n + 1 = 1 − hλ 1 y n ⇒ y n = ( 1 − hλ 1 ) n y 0
So backward Euler is unconditionally stable for decaying problems — pick h h h accuracy , not survival.
For linear f f f nonlinear f f f y n + 1 = y n + h f ( t n + 1 , y n + 1 ) y_{n+1}=y_n+hf(t_{n+1},y_{n+1}) y n + 1 = y n + h f ( t n + 1 , y n + 1 ) Newton's method at each step:
g ( y ) = y − y n − h f ( t n + 1 , y ) = 0 , y ( k + 1 ) = y ( k ) − g ( y ( k ) ) g ′ ( y ( k ) ) , g ′ = 1 − h ∂ f ∂ y . g(y)=y-y_n-h f(t_{n+1},y)=0,\qquad y^{(k+1)}=y^{(k)}-\frac{g(y^{(k)})}{g'(y^{(k)})},\quad g'=1-h\frac{\partial f}{\partial y}. g ( y ) = y − y n − h f ( t n + 1 , y ) = 0 , y ( k + 1 ) = y ( k ) − g ′ ( y ( k ) ) g ( y ( k ) ) , g ′ = 1 − h ∂ y ∂ f .
Intuition WHY Newton, not simple iteration?Simple fixed-point iteration y ( k + 1 ) = y n + h f ( t n + 1 , y ( k ) ) y^{(k+1)}=y_n+hf(t_{n+1},y^{(k)}) y ( k + 1 ) = y n + h f ( t n + 1 , y ( k ) ) ∣ h ∂ f / ∂ y ∣ < 1 |h\partial f/\partial y|<1 ∣ h ∂ f / ∂ y ∣ < 1 stiff restriction we wanted to escape! Newton has no such limit, so we keep big steps.
Worked example Example 1 — forward vs backward,λ = − 100 \lambda=-100 λ = − 100 h = 0.05 h=0.05 h = 0.05 y 0 = 1 y_0=1 y 0 = 1 Forward: 1 + h λ = 1 − 5 = − 4 1+h\lambda = 1-5 = -4 1 + hλ = 1 − 5 = − 4 y n = ( − 4 ) n y_n=(-4)^n y n = ( − 4 ) n 1 , − 4 , 16 , − 64 , … 1,-4,16,-64,\dots 1 , − 4 , 16 , − 64 , … blows up & oscillates .
Why this step? ∣ 1 + h λ ∣ = 4 > 1 |1+h\lambda|=4>1 ∣1 + hλ ∣ = 4 > 1
Backward: 1 1 − h λ = 1 1 + 5 = 1 6 \dfrac{1}{1-h\lambda}=\dfrac{1}{1+5}=\dfrac16 1 − hλ 1 = 1 + 5 1 = 6 1 y n = ( 1 / 6 ) n y_n=(1/6)^n y n = ( 1/6 ) n 1 , 0.167 , 0.028 , ⋯ → 0 1,0.167,0.028,\dots\to0 1 , 0.167 , 0.028 , ⋯ → 0 Why this step? ∣ 1 / ( 1 − h λ ) ∣ = 1 / 6 < 1 |1/(1-h\lambda)|=1/6<1 ∣1/ ( 1 − hλ ) ∣ = 1/6 < 1 h h h e − 100 t → 0 e^{-100t}\to0 e − 100 t → 0
Worked example Example 2 — one backward-Euler step on a nonlinear ODESolve y ′ = − 50 y + sin t y'=-50\,y + \sin t y ′ = − 50 y + sin t y ( 0 ) = 1 y(0)=1 y ( 0 ) = 1 h = 0.1 h=0.1 h = 0.1 y 1 y_1 y 1 y y y y 1 = y 0 + h ( − 50 y 1 + sin t 1 ) y_1 = y_0 + h(-50 y_1 + \sin t_1) y 1 = y 0 + h ( − 50 y 1 + sin t 1 ) t 1 = 0.1 t_1=0.1 t 1 = 0.1 y 1 = 1 + 0.1 ( − 50 y 1 + sin 0.1 ) y_1 = 1 + 0.1(-50 y_1 + \sin 0.1) y 1 = 1 + 0.1 ( − 50 y 1 + sin 0.1 ) y 1 + 5 y 1 = 1 + 0.1 ( 0.0998 ) y_1 + 5y_1 = 1 + 0.1(0.0998) y 1 + 5 y 1 = 1 + 0.1 ( 0.0998 ) 6 y 1 = 1.00998 ⇒ y 1 = 0.16833 6 y_1 = 1.00998 \Rightarrow y_1 = 0.16833 6 y 1 = 1.00998 ⇒ y 1 = 0.16833 Why this step? We moved the − 5 y 1 -5y_1 − 5 y 1 y 1 y_1 y 1 f f f
Worked example Example 3 — choosingh h h y ′ = − 1000 y y'=-1000y y ′ = − 1000 y h < 2 / 1000 = 0.002 h<2/1000=0.002 h < 2/1000 = 0.002 t = 10 t=10 t = 10 ≥ 5000 \geq 5000 ≥ 5000 just
for stability . Backward Euler with h = 0.1 h=0.1 h = 0.1 100 100 100 Why? No stability cap — only accuracy decides.
Common mistake "Implicit = more accurate."Why it feels right: implicit methods are more expensive and handle stiff problems, so they seem superior.
Fix: Backward Euler is still only first-order accurate (O ( h ) O(h) O ( h ) stability , not accuracy. Accuracy and stability are different axes.
Common mistake "Smaller step always means more stable, so just shrinkh h h Why it feels right: for explicit methods shrinking h h h Fix: True, but stiffness makes the required h h h h h h
Common mistake Confusing thesign in 1 − h λ 1-h\lambda 1 − hλ Why it feels right: people copy the forward formula 1 + h λ 1+h\lambda 1 + hλ Fix: Backward Euler evaluates f f f ( 1 − h λ ) y n + 1 = y n (1-h\lambda)y_{n+1}=y_n ( 1 − hλ ) y n + 1 = y n λ < 0 \lambda<0 λ < 0 1 − h λ > 1 1-h\lambda>1 1 − hλ > 1 y y y
Recall Feynman: explain to a 12-year-oldImagine pushing a swing that calms down super fast. If you only look at it and guess where it'll be
a moment later (forward Euler), and you wait too long between looks, your guess overshoots wildly and
you "fall off." Backward Euler is smarter: you ask "where must I be so that things stay calm AT the
next moment?" and solve for that. Because you aim at the calm future, you never overshoot — so you can
take big, lazy steps and still get it right.
"Future-aim stays tame." Backward Euler aims f f f future point t n + 1 t_{n+1} t n + 1 tame (stable).
Sign trick: B ackward → B ig denominator 1 − h λ 1-h\lambda 1 − hλ y y y
What makes an ODE stiff? It has solution components on very different time scales, forcing explicit methods to use tiny
h h h for stability even where the solution is smooth.
Forward Euler formula (test problem y ′ = λ y y'=\lambda y y ′ = λ y y n + 1 = ( 1 + h λ ) y n y_{n+1}=(1+h\lambda)y_n y n + 1 = ( 1 + hλ ) y n .
Forward Euler stability condition? ∣ 1 + h λ ∣ < 1 |1+h\lambda|<1 ∣1 + hλ ∣ < 1 ; for real
λ < 0 \lambda<0 λ < 0 ,
h < 2 / ∣ λ ∣ h<2/|\lambda| h < 2/∣ λ ∣ .
Backward Euler formula? y n + 1 = y n + h f ( t n + 1 , y n + 1 ) y_{n+1}=y_n+h f(t_{n+1},y_{n+1}) y n + 1 = y n + h f ( t n + 1 , y n + 1 ) (implicit).
Backward Euler amplification on y ′ = λ y y'=\lambda y y ′ = λ y y n + 1 = 1 1 − h λ y n y_{n+1}=\dfrac{1}{1-h\lambda}y_n y n + 1 = 1 − hλ 1 y n .
Backward Euler stability condition? ∣ 1 / ( 1 − h λ ) ∣ < 1 |1/(1-h\lambda)|<1 ∣1/ ( 1 − hλ ) ∣ < 1 , true for all
h > 0 h>0 h > 0 when
Re λ < 0 \operatorname{Re}\lambda<0 Re λ < 0 .
What does A-stable mean? The stability region contains the entire left half-plane
Re ( h λ ) < 0 \operatorname{Re}(h\lambda)<0 Re ( hλ ) < 0 .
Order of accuracy of backward Euler? First order,
O ( h ) O(h) O ( h ) — same as forward Euler.
Why use Newton's method inside backward Euler? To solve the nonlinear implicit equation without the step-size restriction that fixed-point iteration would impose.
Implicit vs explicit defining difference? Explicit uses
f f f at known
t n t_n t n ; implicit uses
f f f at unknown
t n + 1 t_{n+1} t n + 1 , so
y n + 1 y_{n+1} y n + 1 must be solved for.
Forward Euler method Region of absolute stability A-stability and L-stability Trapezoidal / Crank–Nicolson method Newton's method for root finding Runge–Kutta methods Eigenvalues and time scales of linear systems
Wildly different time scales
Backward difference at t n+1
abs 1 over 1-h lambda < 1
Intuition Hinglish mein samjhoDekho, kuch ODEs "stiff" hoti hain — matlab unme ek part bahut tezi se decay karta hai (jaise
e − 1000 t e^{-1000t} e − 1000 t h h h stability ke liye, accuracy ke liye nahi. y ′ = λ y y'=\lambda y y ′ = λ y y n + 1 = ( 1 + h λ ) y n y_{n+1}=(1+h\lambda)y_n y n + 1 = ( 1 + hλ ) y n ∣ 1 + h λ ∣ < 1 |1+h\lambda|<1 ∣1 + hλ ∣ < 1 λ = − 1000 \lambda=-1000 λ = − 1000 h < 0.002 h<0.002 h < 0.002
Iska ilaaj hai backward (implicit) Euler : f f f naye point t n + 1 t_{n+1} t n + 1 y n + 1 = y n + h f ( t n + 1 , y n + 1 ) y_{n+1}=y_n+h f(t_{n+1},y_{n+1}) y n + 1 = y n + h f ( t n + 1 , y n + 1 ) y n + 1 y_{n+1} y n + 1 y n + 1 = 1 1 − h λ y n y_{n+1}=\frac{1}{1-h\lambda}y_n y n + 1 = 1 − hλ 1 y n λ < 0 \lambda<0 λ < 0 har h h h A-stable kehte hain — poora left half-plane cover karta hai.
Ek important baat: backward Euler zyada accurate nahi hai, woh bhi sirf first-order (O ( h ) O(h) O ( h ) stability hai — ab tum step size accuracy ke hisaab se choose karo, survival ke
liye nahi. Nonlinear case me implicit equation ko Newton's method se solve karte hain, kyunki simple
iteration phir wahi stiff restriction laga deti hai. Yaad rakho: "Future-aim stays tame" — future point
pe aim karo, system tame (stable) rahega.
