Handout19

Handout19 - Course 18.327 and 1.130 Wavelets and Filter...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Course 18.327 and 1.130 Wavelets and Filter Banks Numerical solution of PDEs: Galerkin approximation; wavelet integrals (projection coefficients, moments and connection coefficients); convergence Numerical Solution of Differential Equations Main idea: look for an approximate solution that lies in Vj. Approximate solution should converge to true solution as j → ∞. Consider the Poisson equation · § leave boundary ∂2µ = f(x) --------------- © conditions till later ¹ ∂x2 Approximate solution: uapprox(x) = ¦ c[k]2j/2 φ(2j x – k) ----------k φj,k(x) trial functions 2 1 Method of weighted residuals: Choose a set of test functions, gn(x), and form a system of equations (one for each n). ³ ∂2uapprox ∂x2 gn(x)dx = ³ f(x)gn(x) dx One possibility: choose test functions to be Dirac delta functions. This is the collocation method. gn(x) = δ(x – n/2j) ¦ c[k]φj,k(n/2j) = f(n/2j) ″ Ÿ n integer ---------------------------- k 3 Second possibility: choose test functions to be scaling functions. • Galerkin method if synthesis functions are used (test functions = trial functions) • Petrov-Galerkin method if analysis functions are used e.g. Petrov-Galerkin ~ gn(x) = φj,n(x) Ÿ ∞ ∂2 2 -∞ ∂x ¦ c[k] ³ k ∈ ~ Vj ∞ ~ ~ φj,k(x) . φj,n(x) dx = ³ f(x)φj,n(x) dx -------- Note: Petrov-Galerkin -∞ ≡ Galerkin in orthogonal case 4 2 Two types of integrals are needed: (a) Connection Coefficients ∞ ∞ ~ 2 ~ ³ ∂∂x2 φj,k(x) . φj,n(x)dx = 22j ³ 2j/2φ″(2jx - k)2j/2φ(2jx - n)dx -∞ -∞ ∞ ~ = 22j ³ φ″(τ)φ(τ + k – n) dτ φ″ -∞ = 22jh∂2/∂x2 [n – k] where h∂2/∂x2 [n] is defined by ∞ ~ h∂2/∂x2 [n] = ³ φ″(t)φ(t – n)dt -----------------∞ connection coefficients 5 (b) Expansion coefficients ∞ ~ The integrals ³ f(x)φj,n(x)dx are the coefficents for -∞ the expansion of f(x) in Vj. fj(x) = ¦ rj[k] φj,k(x) -------------------------- k with ∞ rj[k] = ³ f(x) ~j,k(x) dx φ -∞ -------------------------- So we can write the system of Galerkin equations as a convolution: 22j ¦ c[k]h∂2/∂x2[n – k] = rj[n] ---------------k 6 3 ŸSolve a deconvolution problem to find c[k] and then find uapprox using equation . Note: we must allow for the fact that the solution may be non-unique, i.e. H∂2/∂x2(ω) may have zeros. Familiar example: 3-point finite difference operator h∂2/∂x2[n] = {1, -2, 1} H∂2/∂x2(z) = 1 –2z–1 + z–2 = (1 – z–1)2 Ÿ H∂2/∂x2(ω) has a 2nd order zero at ω = 0. Suppose u0(x) is a solution. Then u0(x) + Ax + B is also a solution. Need boundary conditions to fix uapprox(x). 7 Determination of Connection Coefficients ∞ ~ h∂2/∂x2[n] = ³ φ ″(t) φ(t – n)dt -∞ Simple numerical quadrature will not converge if φ ″ (t) behaves badly. Instead, use the refinement equation to formulate an eigenvalue problem. φ(t) = 2 ¦ f0[k]φ(2t – k) k φ″(t) = 8 ¦ f0[k]φ″(2t – k) k ~ ~ φ(t – n) = 2 ¦ h0[φ(2t – 2n - ) Multiply and Integrate  So h∂2/∂x2[n] = 8 ¦ f0[k] ¦ h0[h∂2/∂x2[2n +  - k] k  8 4 Daubechies 6 scaling function First derivative of Daubechies 6 scaling function 9 Reorganize as h∂2/∂x2[n] = 8¦ h0[m – 2n](¦ f0[m – k]h∂2/∂x2[k]) m k Matrix form h∂2/∂x2 = 8 A B h∂2/∂x2 m = 2n + eigenvalue problem Need a normalization condition use the moments of the scaling function: If h0[n] has at least 3 zeros at π, we can write ∞ ~ ¦ µ2[k]φ(t – k) = t2 ; µ2[k] = ³ t2φ(t – k)dt k -∞ ~ Differentiate twice, multiply by φ(t) and integrate: ¦ µ2[k]h∂2/∂x2[- k] = 2! k Normalizing condition 10 5 Formula for the moments of the scaling function ∞  µ k _ ³τφ(τ - k)dτ -∞ ∞ Recursive formula ∞ µ0 = ³ φ(τ)dτ = 1 0 µr = 0 -∞  µk = r-1 N 1 r r – i)µ i 0 2r - 1 ¦ ( i ) (¦ h0[k]k k=0 i=0  r  ¦ ( r )k-r µ 0 r=0 11 How to enforce boundary conditions? One idea – extrapolate a polyp-1 nomial: u(x) = ¦ c[k]φj,k(x) = ¦ a[x k =0 Relate c[k] to a[ through moments. Extend c[k] by extending underlying polynomial. Extrapolated polynomial should satisfy boundary constraints: Dirichlet: p-1  u(x0) = α Ÿ ¦ a[x0 = α =0 Neumann: p-1 u'(x0) = ⠟ ¦ a[x0-1 = β =0 Constraint on a[ 12 6 Convergence Synthesis scaling function: φ(x) = 2 ¦ f0[k]φ(2x – k) k We used the shifted and scaled versions, φj,k(x), to synthesize the solution. If F0(ω) has p zeros at π, then we can exactly represent solutions which are degree p – 1 polynomials. In general, we hope to achieve an approximate solution that behaves like u(x) = ¦ c[k]φj,k(x) + O(hp) k where h= 1 2j = spacing of scaling functions 13 Reduction in error as a function of h 14 7 Multiscale Representation e.g. ∂2u/∂x2 = f Expand as J u = ¦ ckφ(x – k) + ¦ ¦ dj,kw(2j x-k) j=0 k k Galerkin gives a system Ku = f with typical entries ∞ ∂2 Km,n = 22j ³ ∂∂x2 w(x – n)w(x-m)dx -∞ 15 Effect of Preconditioner • • Multiscale equations: (WKWT)(Wu) = Wf Preconditioned matrix: Kprec = DWKWTD Simple diagonal preconditioner ª1 «1 « « « D=« « « « « ¬ 1 2 1 2 1 4 1 4 º » » » » » » » » 1» 4¼ 16 8 Matlab Example Numerical solution of Partial Differential Equations 17 The Problem 1. Helmholtz equation: uxx + a u = f • p=6; % Order of wavelet scheme (pmin=3) • a=0 • L = 3; % Period. • nmin = 2; % Minimum resolution • nmax = 7; % Maximum resolution 18 9 Solution at Resolution 2 19 Solution at Resolution 3 20 10 Solution at Resolution 4 21 Solution at Resolution 5 22 11 Solution at Resolution 6 23 Solution at Resolution 7 24 12 Convergence Results >> helmholtz slope = 5.9936 25 13 ...
View Full Document

{[ snackBarMessage ]}

Ask a homework question - tutors are online