Machine epsilon
fl(x)=x(1+ϵ) where ∣ϵ∣≤u ∣xfl(x)−x∣=∣ϵ∣≤u is called relative error. Cancellations occur when subtracting nearby number containing roundoff.
Taylor series
f(x) En+1 ∣En+1∣≤chn+1 =k=0∑infk!f(k)(c)(x−c)k=(n+1)!f(n+1)(ξ)(h:=x−c)n+1Polynomial Interpolation
v(x)= Linear system: j=0∑ncjϕj(x) →linearly independent iff v(x)=0 ∀ x→cj=0 ∀ j)ϕ0(x0)ϕ0(x1)⋮ϕ0(xn)ϕ1(x0)ϕ1(x1)⋮ϕ1(xn)⋯⋯⋱⋯ϕn(x0)ϕn(x1)⋮ϕn(xn)c0c1⋮cn=y0y1⋮yn Monomial basis: X: if ϕj(x)=xj, j=0,1,...,n →v(x)=j=0∑ncjxjpn(xi)=c0+c1xi+c2xi2+⋯+cnxin=yiVandermonde matrix→det(X)=i=0∏n−1[j=i+1∏n(xj−xi)]xi are distinct:∙ det(X)=0∙ X is nonsingular∙ system has unique solution∙ unique polynomial of degree≤n that interpolates the data∙ can be poorly conditioned, work is O(n3) Lagrange basis: Lj(xi)={01if i=jif i=jLj(x)=i=0,i=j∏nxj−xix−xipn(xi)=j=0∑nyjLj(xi)=j=0∑i−1yjLj(xi)+yiLi(xi)+j=i+1∑nyjLj(xi)=yi Newton’s basis: ϕj(x)=i=0∏j−1(x−xi),j=0:npn(xi)=c0+c1(xi−x0)+⋯+cn(xi−x0)(xi−x1)⋯(xi−xn−1)=f(xi) Divided differences: f[xi,⋯,xj]=xj−xif[xi+1,⋯,xj]−f[xi,⋯,xj−1]∙ at x=x0 then c0=f(x0)=f[x0]∙ at x=x1 then c1=x1−x0f(x1)−f(x0)=f[x0,x1]∙ at x=x2 then c2=(x2−x0)(x2−x1)f(x2)−c0−c1(x2−x0)=x2−x0x2−x1f(x2)−f(x1)−x1−x0f(x1)−f(x0)=f[x0,x1,x2]∴∀x∈[a,b] ∃ ξ=ξ(x)∈(a,b) : f(x)−pn(x)=(n+1)!fn+1(ξ)i=0∏n(x−xi)∴ Error: ∣f(x)−pn(x)∣≤4(n+1)Mhn+1where: ∙ M=maxa≤t≤b∣fn+1(t)∣∙ h=nb−a∙ xi=a+ih for i=0,1,⋯,nBasic Numeric Integration
If=∫abf(x)dx≈j=0∑najf(xj) (quadrature rule)∙ x0,⋯,xn be distinct points in [a,b]∙ pn(x) be interpolating polynomial of f→ ∫abf(x)dx≈∫abpn(x)dx∙ Uses Lagrange form: ∫abf(x)dx≈j=0∑nf(xj)∫abLj(x)dx=j=0∑nf(xj)aj Trapezoidal rule: ∴ Error: then: From MVT: ∴ f(x)≈p1(x)=f(x0)L0(x)+f(x1)L1(x) (n=1,x0=a,x1=b)If=∫abf(x)dx≈f(a)∫aba−bx−bdx+f(b)∫abb−ax−adx =2b−a[f(a)+f(b)]f(x)−p1(x)=21f′′(ξ(x))(x−a)(x−b)∫ab(f(x)−p1(x))dx=21∫abf′′(ξ(x))(x−a)(x−b)dx∃ η∈(a,b) : ∫abf′′(ξ(x))(x−a)(x−b)dx=f′′(η)∫ab(x−a)(x−b)dxError of Trapezoidal rule: If−Itrap=−12f′′(η)(b−a)3 Midpoint rule: ∴ ∴ If≈Imid=(b−a)f(2a+b)Let m=2a+b→f(x)=f(m)+f′(m)(x−m)+21f′′(ξ(x))(x−m)2If=∫abf(x)=(b−a)f(m)+21∫abf′′(ξ(x))(x−m)2dx∃ η∈(a,b) : 21∫abf′′(ξ(x))(x−m)2dx=24f′′(η)(b−a)3Error of Midpoint rule: If−Imid=24f′′(η)(b−a)3 Simpson’s rule: ∴ If≈Isimp=6b−a[f(a)+4f(2a+b)+f(b)](p2(x),n=2,x0=a,x1=2a+b,x2=b)Error of Simpson’s rule: If−ISimpson=−90f(4)(η)(2b−a)5, η∈(a,b)Composite Numeric Integration
∙ subdivide [a,b] int r subintervals∙ h=rb−a length per interval∙ ti=a+ih for i=0,1,⋯,rt0=a,tr=b → ∫abf(x)dx=i=1∑r∫ti−1tif(x)dx Composite Trapezoidal rule: Error: Composite Simpson rule: Error: Composite Midpoint rule: Error: Icf=2h[f(a)+f(b)]+hi=1∑r−1f(ti)If−Icf=−12f′′(μ)(b−a)h2Ics=3h[f(a)+2i=1∑r/2−1f(t2i)+4i=1∑r/2f(t2i−1)+f(b)]If−Ics=−180f(4)(ζ)(b−a)h4Icm=hi=1∑rf(a+(i−1/2)h)If−Icm=−24f′′(η)(b−a)h2Linear Least Squares
_Find cj such that ∑k=0m(v(x∗k)−yk)2=∑∗k=0m(∑∗j=0ncjϕj(xk)−yk)2 is minimised* Conditions: ∂a∂ϕ=0,∂b∂ϕ=0
Linear fit: yk [∑k=0mxk2∑k=0mxk∑k=0mxkm+1][ab] p →[sppm+1][ab] ↔Az =axk+b,k=1,⋯,m=[∑k=0mxkyk∑k=0myk]=k=0∑mxk,q=k=0∑myk,r=k=0∑mxkyk,s=k=0∑mxk2=[rq]=x0x1⋮xm11⋮1[ab]=y0y1⋮ym=f is overdetermined Solving linear system: r ∣∣r∣∣22 Let ϕ(x) Conditions:∂xk∂ϕ 0 →i=1∑maikj=1∑naijxj =b−Ax=i=1∑mri2=i=1∑m(bi−j=1∑naijxj)2=21∥r∥22=21i=1∑m(bi−j=1∑naijxj)2=0,k=1,⋯,n=i=1∑m(bi−j=1∑naijxj)(−aik)=i=1∑maikbi,k=1,⋯,n (equivalent to ATAx=ATb) ATAx If A has a full-column rank, x =ATb is called the normal equationsxmin∥b−Ax∥2 has uniq sol:=(ATA)−1ATb=A+b Adaptive Simpson: find I S1=S(a,b) E1=E(a,b) Q : ∣Q−I∣≤tol=∫abf(x)dx=S(a,b)+E(a,b)=6h[f(a)+4f(2a+b)+f(b)]=−901(2h)5f(4)(ξ),ξ between a and b S= quadSimpson(f,a,b,tol)h=b−a,c=2a+bS1=6h[f(a)+4f(2a+b)+f(b)]S2=12h[f(a)+4f(2a+c)+2f(c)+4f(2c+b)+f(b)]E~2=151(S2−S1)if∣E~2∣≤tol return Q=S2+E~2else Q1=quadSimpson(f,a,c,tol/2) Q2=quadSimpson(f,c,b,tol/2) return Q=Q1+Q2Newton’s Method for Nonlinear equations
xn+1=xn−f′(xn)f(xn)
Convergence: if f,f′,f′′ are continuous in a neighborhood of a root r of f and f′(r)=0, then ∃δ >0 : ∣r−x0∣≤δ, then ∀xn : :∣r−xn∣≤δ,∣r−xn+1∣≤c(δ)∣r−xn∣2 ∣en+1∣≤c(δ)∣en∣2 (Quadratic convergence, order is 2)
Let c(δ)=21∗min∣r−x∣≤δ∣f′(x)∣max∣r−x∣≤δ∣f′′(x)∣
For linear system: denote x=(x1,x2,⋯,xn)T and F=(f1,f2,⋯,fn), find x∗ such that F(x∗)=0
F(x(k))+F′(x(k))(x(k+1)−x(k)) F′(x(k)) Let s ∴ F′(x(k))s x(k+1) =0is the Jacobian of F at x(k)=x(k+1)−x(k)=−F(x(k))=x(k)+sIVP in ODEs.
Given y′=f(t,y),y(a)=c, find y(t) for t∈[a,b] y′ System of n first-order: y′ ≡y′(t)≡dtdy=f(t,y),f:R×Rn→Rn Forward Euler’s method (explicit): yti+1 where: h h t0 ≈y(ti)+hf(ti,y(ti))=Nb−a,N>1=step size=a,ti=a+ih,i=1,2,⋯,NBackward Euler’s method (implicit): yi+1=yi+hf(ti+1,yi+1)
FE Stability: y′ Exact sol: y(t) FE sol with constant stepsize h: yi+1 To be numerically stable: h BE Stability: y′ ∣yi+1∣ =λy,y(0)=y0=y0eλt=(1+hλ)yi=(1+hλ)i+1y0≤∣λ∣2=λy,y(0)=y0=∣1−hλ∣1∣yi∣≤∣yi∣ ∀ h>0Non-linear, then apply Newton’s methods
Order, Error, Convergence and Stiffness
Local truncation error of FE: Local truncation error of BE: di=hy(ti+1)−y(ti)−f(ti,y(ti))=2hy′′(ηi) (q=1)di=−2hy′′(ξi) (q=1)A method of order q ifq is the lowest positive int such that any smooth exact sol of y(t):maxi∣di∣=O(hq)
Global error: ei Consider u′ =y(ti)−yi,i=0,1,⋯,N=f(t,u),u(ti−1)=yi−1, local error: li=u(ti) Convergence: imaxei=imax∣y(ti)−yi∣→0 as h→0Stiffness is when the stepsize is restricted by stability rather than accuracy
Runge-Kutta Methods
Implicit trapezoidal: y′(t) yi+1 di=O(h2) Explicit trapezoidal: Y yi+1 di=O(h2) Implicit midpoint: yi+1 Explicit midpoint: Y =f(t,y),y(ti)=yi=yi+2h[f(ti,yi)+f(ti+1,yi+1)]=hy(ti+1)−y(ti)−21[f(ti,y(ti))+f(ti+1,y(ti+1))]=yi+hf(ti,yi)=yi+2h[f(ti,yi)+f(ti+1,Y)]=hy(ti+1)−y(ti)−21[f(ti,y(ti))+f(ti+1,y(ti)+hf(ti,y(ti)))]=yi+hf(ti+h/2,(yi+yi+1)/2)=yi+2hf(ti,yi)Classical RK4: based on Simpson’s quadrature rule, O(h4) accuracy
Y1 Y2 Y3 Y4 yi+1 =yi=yi+2hf(ti,Y1)=yi+2hf(ti+2h,Y2)=yi+hf(ti+2h,Y3)=yi+6h[f(ti,Y1)+2f(ti+2h,Y2)+2f(ti+2h,Y3)+f(ti+1,Y4)]