Numerical methods for a nonlinear parabolic problem Jeffrey M. Connors February 20, 2015 The problem is to find u(x, t) : R2 → R solving ut = uxx + f(u), (t > 0), (1) f(u) = u(u − β)(1 − u), 0 < β < 1 2 , (2) u(−∞, t)=1, (3) u(∞, t)=0, (4) u(x, 0) = u0(x), (traveling wave). (5) We shall describe two numerical methods to estimate the solution. We must truncate the domain for computational purposes, so we estimate the boundary conditions as u(−L, t) = 1 and u(L, t) = 0. Here, L $ 1 is an adjustable parameter (that should be set near the top of the code). Let M be the number of interior nodes xj in space, j = 1, 2,…,M, with uniform grid spacing ∆x = 2L/(M + 1); x0 = −L and xM+1 = L. Here, the goal is for u(xj , t) ≈ uj (t). We discretize (1) in space as follows: ∂ ∂t uj (t) = 1 ∆x2 (uj+1 − 2uj + uj−1) + f(uj ), (6) for j = 1,…M; u0 = 1 and uM+1 = 0. We expect that when ∆x is small, |u(xj , t) − uj (t)| ≈ C∆x2 for some C > 0. Consider first using a time-stepping method such as Crank-Nicolson (below). Let ∆t > 0 be the size of the time step, with N total time steps to be taken. Let our approximations now be un j ≈ u(xj , tn) with u0 j = u0(xj ) and define t n+1/2 = t n + ∆t/2. The method is derived by noting that ut(t n+1/2) = u(t n+1) − u(t n) ∆t + O ! ∆t 2″ and then estimating all other terms in (6) at time t n+1/2. Let un+1/2 j = (un+1 j + un j )/2. There are a number of ways to define the method, but it is convenient to use un+1 j − un j ∆t = 1 ∆x2 # un+1/2 j+1 − 2un+1/2 j + un+1/2 j−1 $ + f(un+1/2 j ), (7) for j = 1, 2,…,M. Note that these are only the interior nodes; boundary values get absorbed into the right-hand side of the equations as data. We are not done yet. We need to address the nonlinearity. Two approaches will be taken for this purpose, as outlined below. Extrapolated Crank-Nicolson We replace the nonlinear term f(un+1/2 j ) in (7) with the second-order accurate extrapolation E(f) n+1/2 j = 3 2 f ! un j ” − 1 2 f ! un−1 j ” . 1 The result is un+1 j − un j ∆t = 1 ∆x2 # un+1/2 j+1 − 2un+1/2 j + un+1/2 j−1 $ + E(f) n+1/2 j , which we manipulate into the form un+1 j − ∆t 2∆x2 ! un+1 j+1 − 2un+1 j + un+1 j−1 ” = un j + ∆t 2∆x2 ! un j+1 − 2un j + un j−1 ” + ∆t E(f) n+1/2 j . (8) This is implemented by solving a linear system at each time step of the form (I + γD)$un+1 = F n = (I − γD)$un + ∆tE$ (f) n+1/2 + 2γ[1 0 … 0]T with I the M × M identity matrix, D a diffusion matrix with γ = ∆t/(2∆x2), $uk a column vector with components uk j and E$ (f)n+1/2 a column vector with components E(f) n+1/2 j . The last term on the right-hand side is due to boundary conditions. Newton’s Method In this case, we iterate to resolve the nonlinearity using Newton’s method to find Φ($un+1) = 0, where Φ($z)=(I + γD)$z − (I − γD)$un − ∆tG(($z + $un)/2) − 2γ[1 0 … 0]T . (9) Here $z has components zj , j = 1,…M and G($η)=[f(η1) … f(ηM)]T . The iterates will be $z(k) , k = 0, 1,…. We will choose the initial guess $z(0) = $un at each time step. Then apply the usual Newton iteration: DΦ($z(k) )∆$z(k) = −Φ($z(k) ) $z(k+1) = $z(k) + ∆$z(k) . One can use (2) and some standard calculus to show that DΦ($z(k) )∆$z(k) = # I + γD − ∆tG($z(k) ) $ ∆$z(k) , G($z(k) ) = diag (g1, g2, …, gM), gj = 1 8 % −3 # z (k) j + un j $2 + 4(1 + β) # z (k) j + un j $ − 4β & , j = 1,…,M. Each iteration requires one to solve a linear system. It is preferable to eliminate ∆$z(k) and solve directly for $z(k+1). This is accomplished via the following equation: # I + γD − ∆tG($z(k) ) $ $z(k+1) =(I − γD)$un + ∆tG(($z(k) + $un)/2) − ∆tG($z(k) )$z(k) + 2γ[1 0 … 0]T . (10) One iterates until ‘$z(k+1) − $z(k) ‘ < &, where & > 0 is a user-defined parameter in the code. In rare instances, it may be necessary to check (in addition) that ‘Φ($z(k+1))’ < &. Finally, Newton’s method may not converge in general, so a maximum k-value, say KMAX, should be specified such that the code exits upon reaching k > KMAX. 2
The post describe two numerical methods to estimate the solution appeared first on Dissertation Help Service.