Review
Download
1 / 36

Review - PowerPoint PPT Presentation


  • 75 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


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