1 / 36

# Review - PowerPoint PPT Presentation

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

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

## PowerPoint Slideshow about ' Review' - alana-howell

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

Paul Heckbert

Computer Science Department

Carnegie Mellon University

15-859B - Introduction to Scientific Computing

15-859B - Introduction to Scientific Computing

y’=-y, stable but slow solution

15-859B - Introduction to Scientific Computing

y’=-y, unstable solution

15-859B - Introduction to Scientific Computing

• 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

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

• 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

• 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

• 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

15-859B - Introduction to Scientific Computing

• 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

• 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

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

• 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

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

15-859B - Introduction to Scientific Computing

dx=.05

dt=.06

Euler’s method unstable when step too large

15-859B - Introduction to Scientific Computing

• 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

• 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

• 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

• 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

• 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

• Generally for symmetric positive definite, only.

• Convert linear system Ax=b

• into optimization problem: minimize xTAx-xTb

• a parabolic bowl

• 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

• 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

• 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

• relax on a given grid a few times

• coarsen (restrict) a grid

• refine (interpolate) a grid

15-859B - Introduction to Scientific Computing

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

• Gauss-Seidel

• converges for all symmetric positive definite A

• convergence rate determined by condition number

• note that condition number typically larger for finer grids

• 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

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)

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

multigrid O(k2) = O(n)

optimal!

15-859B - Introduction to Scientific Computing

• 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

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

• 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

• generality – approximate anything well

• discontinuities, nonperiodicity, ...

• 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

• 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

• 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