Boundary Value Problems, Lecture notes of Numerical Methods in Engineering

Download Boundary Value Problems and more Lecture notes Numerical Methods in Engineering in PDF only on Docsity!

Boundary Value Problems Numerical Methods for BVPs

Outline

(^1) Boundary Value Problems (^2) Numerical Methods for BVPs Michael T. Heath Scientific Computing 2 / 45

Boundary Value Problems Numerical Methods for BVPs Boundary Values Existence and Uniqueness Conditioning and Stability

Boundary Value Problems

Side conditions prescribing solution or derivative values at specified points are required to make solution of ODE unique For initial value problem, all side conditions are specified at single point, say t 0 For boundary value problem (BVP), side conditions are specified at more than one point kth order ODE, or equivalent first-order system, requires k side conditions For ODEs, side conditions are typically specified at endpoints of interval [a, b], so we have two-point boundary value problem with boundary conditions (BC) at a and b. Michael T. Heath Scientific Computing 3 / 45

• Weighted residual technique (WRT):
  • Approximate the solution by an expansion, u(x) := X^ n j= uj j (x), and solve for coecients uj such that the ODE is satisfied in some weighted sense.
  • That is, rather than enforcing r(x) = 0 at isolated points, we require r(x) to be orthogonal to a set of weight functions, (^) i(x): Z (^) b a i(x)^ r(x)^ dx^ = Z (^) b a i(x)^ L(u)^ ^ L(˜u)^ dx^ = 0,^ or Z (^) b a i(x)^ L(u)^ = Z (^) b a i(x)^ L(˜u)^ dx for i = 0, 1, ...
  • Note that if (^) i(x) = (x xi) (Dirac delta function), we recover collocation.
  • Most often, the test-space and trial space are the same: (^) i := i. (Galerkin case.)
  • Finite element, spectral, spectral element methods are examples of WRTs.
  • WRTs have many advantages over collocation in terms of flexibility of basis func- tions, application of boundary conditions, etc., and are generally preferred over collocation.

Boundary Value Problems Numerical Methods for BVPs Shooting Method Finite Difference Method Collocation Method Galerkin Method

Finite Difference Method

Finite difference method converts BVP into system of algebraic equations by replacing all derivatives with finite difference approximations For example, to solve two-point BVP u 00 = f (t, u, u 0 ), a < t < b with BC u(a) = ↵, u(b) = we introduce mesh points ti = a + ih, i = 0, 1 ,... , n + 1, where h = (b a)/(n + 1) We already have y 0 = u(a) = ↵ and yn+1 = u(b) = from BC, and we seek approximate solution value yi ⇡ u(ti) at each interior mesh point ti, i = 1,... , n Michael T. Heath Scientific Computing 20 / 45

Boundary Value Problems Numerical Methods for BVPs Shooting Method Finite Difference Method Collocation Method Galerkin Method

Finite Difference Method, continued

We replace derivatives by finite difference approximations such as u 0 (ti) ⇡ yi+1 yi 1 2 h u 00 (ti) ⇡ yi+1 2 yi + yi 1 h 2 This yields system of equations yi+1 2 yi + yi 1 h 2 = f

ti, yi, yi+1 yi 1 2 h

to be solved for unknowns yi, i = 1,... , n System of equations may be linear or nonlinear, depending on whether f is linear or nonlinear Michael T. Heath Scientific Computing 21 / 45 Error is O(h^2 ) Error is O(h^2 )

Boundary Value Problems Numerical Methods for BVPs Shooting Method Finite Difference Method Collocation Method Galerkin Method

Finite Difference Method, continued

For these particular finite difference formulas, system to be solved is tridiagonal, which saves on both work and storage compared to general system of equations This is generally true of finite difference methods: they yield sparse systems because each equation involves few variables Michael T. Heath Scientific Computing 22 / 45

Example: Convection-Diffusion Equation

Comments About Computing Error Norms

• Be careful with the l 2 vector norm!
• Even though max |ei| ! 0 with n ! 1, we can still have ||e|| grow with n. Why?
• When solving di↵erential equations, it is better to use norms that approximate their continuous counterparts. Thus ||e|| 2 = Z ⌦ e 2 dx 1 2 ⇡ " 1 n X^ n i= |ei| 2

1

2 ||e|| 1 = max ⌦ |e| ⇡ max i |ei|

• The issue can also be resolved by measuring relative error: error := ||e|| ||u|| for some appropriate vector norm.
• Still, best to start with a norm that doesn’t scale with n.

Convection-Diffusion Equation

• The issue is that ⌫ ! 0 is a singular perturbation.
• This is true whenever the highest-order derivative is multiplied by a small constant.
• As the constant goes to zero, the number of boundary conditions changes.
• Here,
  • We go from one boundary condition when ⌫ = 0,
  • to two boundary conditions when ⌫ > 0 (even for ⌫ ⌧ 1).
• An example that is not a singular perturbation is uxx + ✏ ux = 1, u(0) = u(1) = 0, ✏ ! 0. This is called a regular perturbation.

Regular / Singular Perturbations You’re Familiar With

• Another example:
• Consider solutions to the quadratic equation: ax 2 + bx + c = 0. Example 1: x 2

• ✏x = 1 : Two roots as ✏ ! 0. Regular perturbation. Example 2: ✏x 2
• x = 1 : x = 1 2 ✏ ± 1 2 ✏ p 1 + 4✏ x 1 = 1 2 ✏ ⇣p 1 + 4✏ 1 ⌘ = 1 2 ✏ 1 + 2✏ 1 + O(✏ 2 ) = 1 + O(✏). x 2 = 1 2 ✏ (2 + O(✏)) ! 1. Singular perturbation.

Convection-Diffusion Equation

❑ What happens when cL/º À 1 in our numerical example?

Nonlinear Example: The Bratu Equation

• Consider 1D di↵usion with nonlinear right-hand side: d 2 u dx 2 = q(x, u) = e u , u(0) = u(1) = 0.
• Discretizing with finite di↵erences (say), Au = e u .
• Nonlinear system: f (u) = 0 , f (u) = Au e u .

f k i = j= Aij u k j ^ e ui ! (Jk) ij = i @uj = Aij e ui .

• If b = 1 and aj = 2 h 2 e uj , then
• At each iteration, modify the tridiagonal matrix A such that Jk = A + e u k i (^) ij , and solve this tridiagonal system in ⇡ 8 n operations.

fj =

uj 1 2 uj uj+

h

2

+ e

uj

= 0, j = 1,... , n

need the Jacobian (2), which is given by the tridiagonal matrix

J =

h

2

0 B B B B B B B B B @

a 1 b

b a 2 b

b

. b

b an

1 C C C C C C C C C A

aj = 2 + h

2

e

uj

. Note that, as is often the case with systems arising

ations, J is sparse. That is, it has a fixed number of nonzeros per row, indepen

as O(n) nonzeros. Moreover, because this system is tridiagonal, the factor co

bratu1a.m / bratu_lin.m