review
Download
Skip this Video
Download Presentation
Review

Loading in 2 Seconds...

play fullscreen
1 / 36

Review - PowerPoint PPT Presentation


  • 76 Views
  • Uploaded on

Review. Paul Heckbert Computer Science Department Carnegie Mellon University. y’=y. y’=-y, stable but slow solution. y’=-y, unstable solution. Jacobian of ODE. ODE: y’=f(t,y) , where y is n -dimensional Jacobian of f is a square matrix

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about ' Review' - alana-howell


An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
review

Review

Paul Heckbert

Computer Science Department

Carnegie Mellon University

15-859B - Introduction to Scientific Computing

slide2
y’=y

15-859B - Introduction to Scientific Computing

y y stable but slow solution
y’=-y, stable but slow solution

15-859B - Introduction to Scientific Computing

y y unstable solution
y’=-y, unstable solution

15-859B - Introduction to Scientific Computing

jacobian of ode
Jacobian of ODE
  • ODE: y’=f(t,y), where y is n-dimensional
  • Jacobian of f is a square matrix
  • if ODE homogeneous and linear then J is constant and y’=Jy
  • but in general J varies with t and y

15-859B - Introduction to Scientific Computing

stability of ode depends on jacobian
Stability of ODE depends on Jacobian

At a given (t,y) find J(t,y) and its eigenvalues

find rmax = maxi { Re[i(J)] }

if rmax<0, ODE stable, locally

rmax =0, ODE neutrally stable, locally

rmax >0, ODE unstable, locally

15-859B - Introduction to Scientific Computing

stability of numerical solution
Stability of Numerical Solution
  • Stability of numerical solution is related to, but not the same as stability of ODE!
  • Amplification factor of a numerical solution is the factor by which global error grows or shrinks each iteration.

15-859B - Introduction to Scientific Computing

stability of euler s method
Stability of Euler’s Method
  • Amplification factor of Euler’s method is I+hJ
  • Note that it depends on h and, in general, on t & y.
  • Stability of Euler’s method is determined by eigenvalues of I+hJ
  • spectral radius (I+hJ)= maxi | i(I+hJ) |
  • if (I+hJ)<1 then Euler’s method stable
    • if all eigenvalues of hJ lie inside unit circle centered at –1, E.M. is stable
    • scalar case: 0<|hJ|<2 iff stable, so choose h < 2/|J|
  • What if one eigenvalue of J is much larger than the others?

15-859B - Introduction to Scientific Computing

implicit methods
Implicit Methods
  • use information from future time tk+1 to take a step from tk
  • Euler method: yk+1 = yk+f(tk,yk)hk
  • backward Euler method: yk+1 = yk+f(tk+1,yk+1)hk
  • example:
  • y’=Ay f(t,y)=Ay
  • yk+1 = yk+Ayk+1hk
  • (I-hkA)yk+1=yk -- solve this system each iteration

15-859B - Introduction to Scientific Computing

stability of backward euler s method
Stability of Backward Euler’s Method
  • Amplification factor of B.E.M. is (I-hJ)-1
  • B.E.M. is stable independent of h (unconditionally stable) as long as rmax<0, i.e. as long as ODE is stable
  • Implicit methods such as this permit bigger steps to be taken (larger h)

15-859B - Introduction to Scientific Computing

large steps ok with backward euler s method
Large steps OK with Backward Euler’s method

15-859B - Introduction to Scientific Computing

popular ivp solution methods
Popular IVP Solution Methods
  • Euler’s method, 1st order
  • backward Euler’s method, 1st order
  • trapezoid method (a.k.a. 2nd order Adams-Moulton)
  • 4th order Runge-Kutta
  • If a method is pth order accurate then its global error is O(hp)

15-859B - Introduction to Scientific Computing

solving boundary value problems
Solving Boundary Value Problems
  • Shooting Method
    • transform BVP into IVP and root-finding
  • Finite Difference Method
    • transform BVP into system of equations (linear?)
  • Finite Element Method
    • choice of basis: linear, quadratic, ...
    • collocation method
      • residual is zero at selected points
    • Galerkin method
      • residual function is orthogonal to each basis function

15-859B - Introduction to Scientific Computing

2 nd order pde classification
2nd Order PDE Classification
  • We classify conic curve ax2+bxy+cy2+dx+ey+f=0 as ellipse/parabola/hyperbola according to sign of discriminant b2-4ac.
  • Similarly we classify 2nd order PDE auxx+buxy+cuyy+dux+euy+fu+g=0:

b2-4ac < 0 – elliptic (e.g. equilibrium)

b2-4ac = 0 – parabolic (e.g. diffusion)

b2-4ac > 0 – hyperbolic (e.g. wave motion)

For general PDE’s, class can change from point to point

15-859B - Introduction to Scientific Computing

example wave equation
Example: Wave Equation
  • utt = c uxx for 0x1, t0
  • initial cond.: u(0,x)=sin(x)+x+2, ut(0,x)=4sin(2x)
  • boundary cond.: u(t,0) = 2, u(t,1)=3
  • c=1
  • unknown: u(t,x)
  • simulated using Euler’s method in t
  • discretize unknown function:

15-859B - Introduction to Scientific Computing

wave equation numerical solution

k+1

k

k-1

j+1

j

j-1

Wave Equation: Numerical Solution

u0 = ...

u1 = ...

for t = 2*dt:dt:endt

u2(2:n) = 2*u1(2:n)-u0(2:n)

+c*(dt/dx)^2*(u1(3:n+1)-2*u1(2:n)+u1(1:n-1));

u0 = u1;

u1 = u2;

end

15-859B - Introduction to Scientific Computing

wave equation results
Wave Equation Results

15-859B - Introduction to Scientific Computing

poor results when dt too big
Poor results when dt too big

dx=.05

dt=.06

Euler’s method unstable when step too large

15-859B - Introduction to Scientific Computing

pde solution methods
PDE Solution Methods
  • Discretize in space, transform into system of IVP’s
  • Discretize in space and time, finite difference method.
  • Discretize in space and time, finite element method.
    • Latter methods yield sparse systems.
  • Sometimes the geometry and boundary conditions are simple (e.g. rectangular grid);
  • Sometimes they’re not (need mesh of triangles).

15-859B - Introduction to Scientific Computing

pde s and sparse systems
PDE’s and Sparse Systems
  • A system of equations is sparse when there are few nonzero coefficients, e.g. O(n) nonzeros in an nxn matrix.
  • Partial Differential Equations generally yield sparse systems of equations.
  • Integral Equations generally yield dense (non-sparse) systems of equations.
  • Sparse systems come from other sources besides PDE’s.

15-859B - Introduction to Scientific Computing

example pde yielding sparse system
Example: PDE Yielding Sparse System
  • Laplace’s Equation in 2-D: 2u = uxx +uyy = 0
    • domain is unit square [0,1]2
    • value of function u(x,y) specified on boundary
    • solve for u(x,y) in interior

15-859B - Introduction to Scientific Computing

sparse matrix storage
Sparse Matrix Storage
  • Brute force: store nxn array, O(n2) memory
    • but most of that is zeros – wasted space (and time)!
  • Better: use data structure that stores only the nonzeros

col 1 2 3 4 5 6 7 8 9 10...

val 0 1 0 0 1 -4 1 0 0 1...

16 bit integer indices: 2, 5, 6, 7,10

32 bit floats: 1, 1,-4, 1, 1

  • Memory requirements, if kn nonzeros:
    • brute force: 4n2 bytes, sparse data struc: 6kn bytes

15-859B - Introduction to Scientific Computing

an iterative method gauss seidel
An Iterative Method: Gauss-Seidel
  • System of equations Ax=b
  • Solve ith equation for xi:
  • Pseudocode:

until x stops changing

for i = 1 to n

x[i]  (b[i]-sum{ji}{a[i,j]*x[j]})/a[i,i]

  • modified x values are fed back in immediately
  • converges if A is symmetric positive definite

15-859B - Introduction to Scientific Computing

conjugate gradient method
Conjugate Gradient Method
  • Generally for symmetric positive definite, only.
  • Convert linear system Ax=b
  • into optimization problem: minimize xTAx-xTb
    • a parabolic bowl
  • Like gradient descent
    • but search in conjugate directions
    • not in gradient direction, to avoid zigzag problem
  • Big help when bowl is elongated (cond(A) large)

15-859B - Introduction to Scientific Computing

convergence of conjugate gradient method
Convergence ofConjugate Gradient Method
  • If K = cond(A) = max(A)/ min(A)
  • then conjugate gradient method converges linearly with coefficient (sqrt(K)-1)/(sqrt(K)+1) worst case.
  • often does much better: without roundoff error, if A has m distinct eigenvalues, converges in m iterations!

15-859B - Introduction to Scientific Computing

example poisson s equation
Example: Poisson’s Equation
  • Sweep of Gauss-Seidel “relaxes” each grid value to be the average of its four neighbors plus an f offset
  • Many relaxations required to solve this on a fixed grid.
  • Multigrid solves it on a hierarchy of grids.

15-859B - Introduction to Scientific Computing

elements of multigrid method
Elements of Multigrid Method
  • relax on a given grid a few times
  • coarsen (restrict) a grid
  • refine (interpolate) a grid

15-859B - Introduction to Scientific Computing

a common multigrid schedule

final solution

finest grid kk

k/2k/2

...

coarsest grid 11

time

=relax

=interpolate

=restrict

A Common Multigrid Schedule

Full Multigrid V Cycle:

15-859B - Introduction to Scientific Computing

some iterative methods
Some Iterative Methods
  • Gauss-Seidel
    • converges for all symmetric positive definite A
  • Conjugate Gradient (CG) Method
    • convergence rate determined by condition number
    • note that condition number typically larger for finer grids
  • Preconditioned Conjugate Gradient
    • instead of solving Av=f, solve M-1Av=M-1f where M-1 is cheap and M is close to A
    • often much faster than CG, but conditioner M is problem-dependent
  • Multigrid
    • convergence rate is independent of condition number, problem size
    • but algorithm must be tuned for a given problem; not as general as others

note: don’t need matrix A in memory – can compute it on the fly!

15-859B - Introduction to Scientific Computing

cost comparison
Cost Comparison

on 2-D Poisson Equation, kk grid, n=k2 unknowns

METHOD COST

Gaussian Elimination O(k6) = O(n3)

Gauss-Seidel O(k4logk) = O(n2logn)

Conjugate Gradient O(k3) = O(n1.5)

FFT/cyclic reduction O(k2logk) = O(nlogn)

multigrid O(k2) = O(n)

optimal!

15-859B - Introduction to Scientific Computing

function representations
Function Representations
  • sequence of samples (time domain)
    • finite difference method
  • pyramid (hierarchical)
  • polynomial
  • sinusoids of various frequency (frequency domain)
    • Fourier series
  • piecewise polynomials (finite support)
    • finite element method, splines
  • wavelet (hierarchical, finite support)
    • (time/frequency domain)

15-859B - Introduction to Scientific Computing

what are wavelets
What Are Wavelets?

In general, a family of representations using:

  • hierarchical (nested) basis functions
  • finite (“compact”) support
  • basis functions often orthogonal
  • fast transforms, often linear-time

15-859B - Introduction to Scientific Computing

nested function spaces for haar basis
Nested Function Spaces for Haar Basis
  • Let Vj denote the space of all piecewise-constant functions represented over 2j equal sub-intervals of [0,1]
  • Vj has basis

15-859B - Introduction to Scientific Computing

function representations desirable properties
Function Representations –Desirable Properties
  • generality – approximate anything well
    • discontinuities, nonperiodicity, ...
  • adaptable to application
    • audio, pictures, flow field, terrain data, ...
  • compact – approximate function with few coefficients
    • facilitates compression, storage, transmission
  • fast to compute with
    • differential/integral operators are sparse in this basis
    • Convert n-sample function to representation in O(nlogn) or O(n) time

15-859B - Introduction to Scientific Computing

time frequency analysis
Time-Frequency Analysis
  • For many applications, you want to analyze a function in both time and frequency
  • Analogous to a musical score
  • Fourier transforms give you frequency information, smearing time.
  • Samples of a function give you temporal information, smearing frequency.
  • Note: substitute “space” for “time” for pictures.

15-859B - Introduction to Scientific Computing

comparison to fourier analysis
Comparison to Fourier Analysis
  • Fourier analysis
    • Basis is global
    • Sinusoids with frequencies in arithmetic progression
  • Short-time Fourier Transform (& Gabor filters)
    • Basis is local
    • Sinusoid times Gaussian
    • Fixed-width Gaussian “window”
  • Wavelet
    • Basis is local
    • Frequencies in geometric progression
    • Basis has constant shape independent of scale

15-859B - Introduction to Scientific Computing

ad