Chapter 26 The Implicit Euler Can Be Written for This Problem as

CHAPTER 26 26.1 (a) h < 2/100,000 = 2×10−5. (b) The implicit Euler can be written for this problem as ( ) y i +1 = yi + − 100,000 y i +1 + 100 ,000e − xi + 1 − e − xi + 1 h which can be solved for y i +1 = y i + 100,000 e − x i + 1 h − e − xi + 1 h 1 + 100,000h The results of applying this formula for the first few steps are shown below. A plot of the entire solution is also displayed x 0 0.1 0.2 0.3 0.4 0.5 y 0 1.904638 1.818731 1.740819 1.67032 1.606531 2 1 0 0 1 2 26.2 The implicit Euler can be written for this problem as y i +1 = yi + ( 30(sin t i +1 − y i +1 ) + 3 cos t i +1 ) h which can be solved for y i +1 = y i + 30 sin t i +1h + 3 cos t i +1h 1 + 30h The results of applying this formula are tabulated and graphed below. x y x 0 0 0.4 0.444484 0.8 0.760677 y x 1.2 0.952306 1.6 0.993242 2 0.877341 y x 2.4 0.622925 2.8 0.270163 3.2 -0.12525 y 3.6 -0.50089 4 -0.79745 1 0.5 0 0 1 2 3 4 -0.5 -1 26.3 (a) The explicit Euler can be written for this problem as x1,i +1 = x1,i + ( 999 x1,i + 1999 x 2 ,i ) h x 2,i +1 = x 2 ,i + ( − 1000 x 1,i − 2000 x 2 ,i ) h Because the step-size is much too large for the stability requirements, the solution is unstable, t 0 0.05 0.1 0.15 0.2 x1 x2 1 150.9 -7204.2 353186 -1.7E+07 dx1 1 -149 7206 -353184 17305943 dx2 2998 -147102 7207803 -3.5E+08 1.73E+10 -3000 147100 -7207805 3.53E+08 -1.7E+10 (b) The implicit Euler can be written for this problem as x1,i +1 = x 1,i + ( 999 x1,i +1 + 1999 x 2,i +1 ) h x 2,i +1 = x 2 ,i + ( − 1000 x 1,i +1 − 2000 x 2 ,i +1 ) h or collecting terms (1 − 999 h) x1,i +1 − 1999 hx 2 ,i +1 = x1,i 1000hx1,i +1 + (1 + 2000h ) x 2,i + 1 = x 2,i or substituting h = 0.05 and expressing in matrix format  − 48.95 − 99.95  x1,i +1  =  x1,i  101   x 2 ,i +1   x 2 ,i   50     Thus, to solve for the first time step, we substitute the initial conditions for the right-hand side and solve the 2x2 system of equations. The best way to do this is with LU decomposition since we will have to solve the system repeatedly. For the present case, because its easier to display, we will use the matrix inverse to obtain the solution. Thus, if the matrix is inverted, the solution for the first step amounts to the matrix multiplication, {} { }  x1,i +1   1886088 . 186648 . 3.752568 1 x  =  − 0.93371 − 0.9141 1 = − 1.84781    2 ,i +1  For the second step (from x = 0.05 to 0.1), { } {  x1,i +1   1886088 . 186648 .  3.752568 = 3.62878 x  =  − 0.93371 − 0.9141 − 184781 . − 181472 .   2 , i + 1   } The remaining steps can be implemented in a similar fashion to give t 0 0.05 0.1 0.15 0.2 x1 1 3.752568 3.62878 3.457057 3.292457 x2 1 -1.84781 -1.81472 -1.72938 -1.64705 The results are plotted below, along with a solution with the explicit Euler using a step of 0.0005. x1 4 2 0 0 0.1 -2 0.2 x2 26.4 First step: Predictor: y10=5.222138+[−0.5(4.143883)+e−2]1 = 3.285532 Corrector: 1 y 1 = 4.143883 + − 0.5( 4.143883) + e −2 − 0.5( 3.285532 ) + e −2.5 0.5 = 3.269562 2 The corrector can be iterated to yield j yi+1j ε a ,% 1 2 3.269562 3.271558 0.061 Second step: Predictor: y20=4.143883+[−0.5(3.271558)+e−2.5]1 = 2.590189 Predictor Modifier: y20 = 2.590189+4/5(3.271558-3.285532) = 2.579010 Corrector: y 21 = 3.271558 + − 0.5( 3.271558 ) + e −2.5 − 0.5(2.579010) + e −3 0.5 = 2.573205 2 The corrector can be iterated to yield 26.5 j yi+1j ε a ,% 1 2 2.573205 2.573931 0.0282 predictor = 3.270674927 Corrector Iteration x y ea 2.5 3.274330476 1.12E-01 2.5 3.273987768 1.05E-02 2.5 3.274019897 9.81E-04 predictor = 2.576436209 Corrector Iteration x y 3 2.57830404 3 2.578128931 3 2.377128276 3 2.377202366 (b) predictor = 0.669232229 Corrector Iteration ea 7.24E-02 6.79E-03 3.32E-02 3.12E-03 x 4.5 4.5 4.5 y 0.666462335 0.666577747 0.666572938 ea 4.16E-01 1.73E-02 7.21E-04 et 0.030654791 0.013342954 0.014064281 predictor = 0.601036948 Corrector Iteration x y ea 5 0.599829531 2.01E-01 5 0.599874809 7.55E-03 26.8 predictor Corrector x 2 2 2 2 predictor Corrector x = 0.737731653 Iteration y 0.660789924 0.665598782 0.665298229 0.665317013 = 0.585786168 Iteration y et 0.028411529 0.020865171 ea 1.16E+01 7.22E-01 4.52E-02 0.002823406 ea 2.5 2.5 2.5 26.9 0.569067395 0.569963043 0.569915062 2.94E+00 1.57E-01 8.42E-03 Option Explicit Sub Dim Dim Dim SimpImplTest() i As Integer, m As Integer xi As Single, yi As Single, xf As Single, dx As Single, xout As Single xp(200) As Single, yp(200) As Single 'Assign values yi = 0 xi = 0 xf = 0.4 dx = 0.05 xout = 0.05 'Perform numerical Integration of ODE Call ODESolver(xi, yi, xf, dx, xout, xp(), yp(), m) 'Display results Sheets("Sheet1").Select Range("a5:b205").ClearContents Range("a5").Select For i = 0 To m ActiveCell.Value = xp(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = yp(i) ActiveCell.Offset(1, -1).Select Next i Range("a5").Select End Sub Sub ODESolver(xi, yi, xf, dx, xout, xp, yp, m) 'Generate an array that holds the solution Dim x As Single, y As Single, xend As Single Dim h As Single m = 0 xp(m) = xi yp(m) = yi x = xi y = yi Do 'Print loop xend = x + xout If (xend > xf) Then xend = xf 'Trim step if increment exceeds end h = dx Call Integrator(x, y, h, xend) m = m + 1 xp(m) = x yp(m) = y If (x >= xf) Then Exit Do Loop End Sub Sub Integrator(x, y, h, xend) Dim ynew As Single Do 'Calculation loop If (xend - x < h) Then h = xend - x Call SimpImpl(x, y, h, ynew) y = ynew If (x >= xend) Then Exit Do Loop End Sub 'Trim step if increment exceeds end Sub SimpImpl(x, y, h, ynew) 'Implement Simple Implicit method ynew = (y + h * FF(x)) / (1 + 1000 * h) x = x + h End Sub Function FF(t) 'Define Forcing Function FF = 3000 - 2000 * Exp(-t) End Function 26.10 All linear systems are of the form dy1 = a11 y1 + a12 y 2 + F1 dt dy2 = a 21 y1 + a22 y2 + F2 dt As shown in the book (p. 730), the implicit approach amounts to solving − a12   y1, i +1   y1,i + F1h  1 − a11 h  − a 21 1 − a22 h   y 2 ,i +1  =  y 2 ,i + F2 h      Therefore, for Eq. 26.6: a11 = −5, a12 = 3, a21 = 100, a22 = −301, F1 =, and F2 = 0, − 3   y1, i +1  =  y1, i  1 + 5h − 100 1 + 301h   y2, i +1   y2 , i       A VBA program written in these terms is Option Explicit Sub Dim Dim Dim StiffSysTest() i As Integer, m As Integer, n As Integer, j As Integer xi As Single, yi(10) As Single, xf As Single, dx As Single, xout As Single xp(200) As Single, yp(200, 10) As Single 'Assign values n = 2 xi = 0 xf = 0.4 yi(1) = 52.29 yi(2) = 83.82 dx = 0.05 xout = 0.05 'Perform numerical Integration of ODE Call ODESolver(xi, yi(), xf, dx, xout, xp(), yp(), m, n) 'Display results Sheets("Sheet1").Select Range("a5:n205").ClearContents Range("a5").Select For i = 0 To m ActiveCell.Value = xp(i) For j = 1 To n ActiveCell.Offset(0, 1).Select ActiveCell.Value = yp(i, j) Next j ActiveCell.Offset(1, -n).Select Next i Range("a5").Select End Sub Sub ODESolver(xi, yi, xf, dx, xout, xp, yp, m, n) 'Generate an array that holds the solution Dim i As Integer Dim x As Single, y(10) As Single, xend As Single Dim h As Single m = 0 x = xi 'set initial conditions For i = 1 To n y(i) = yi(i) Next i 'save output values xp(m) = x For i = 1 To n yp(m, i) = y(i) Next i Do 'Print loop xend = x + xout If (xend > xf) Then xend = xf 'Trim step if increment exceeds end h = dx Call Integrator(x, y(), h, n, xend) m = m + 1 'save output values xp(m) = x For i = 1 To n yp(m, i) = y(i) Next i If (x >= xf) Then Exit Do Loop End Sub Sub Integrator(x, y, h, n, xend) Dim j As Integer Dim ynew(10) As Single Do 'Calculation loop If (xend - x < h) Then h = xend - x Call StiffSys(x, y, h, n, ynew()) For j = 1 To n y(j) = ynew(j) Next j If (x >= xend) Then Exit Do Loop End Sub 'Trim step if increment exceeds end Sub StiffSys(x, y, h, n, ynew) Dim j As Integer Dim FF(2) As Single, b(2, 2) As Single, c(2) As Single, den As Single Call Force(x, FF()) 'MsgBox "pause" b(1, 1) = 1 + 5 * h b(1, 2) = -3 * h b(2, 1) = -100 * h b(2, 2) = 1 + 301 * h For j = 1 To n c(j) = y(j) + FF(j) * h Next j den = b(1, 1) * b(2, 2) - b(1, 2) * b(2, 1) ynew(1) = (c(1) * b(2, 2) - c(2) * b(1, 2)) / den ynew(2) = (c(2) * b(1, 1) - c(1) * b(2, 1)) / den x = x + h End Sub Sub Force(t, FF) 'Define Forcing Function FF(0) = 0 FF(1) = 0 End Sub The result compares well with the analytical solution. If a smaller step size were used, the solution would improve 26.11 (Errata for first printing) Last sentence of problem statement should read: Test the program by duplicating Example 26.4. Later printings should have this correction. Option Explicit Sub Dim Dim Dim Dim Dim NonSelfStartHeun() n As Integer, m As Integer, i As Integer, iter As Integer iterp(1000) As Integer xi As Single, xf As Single, yi As Single, h As Single x As Single, y As Single xp(1000) As Single, yp(1000) As Single xi = -1 xf = 4 yi = -0.3929953 n = 5 h = (xf - xi) / n x = xi y = yi m = 0 xp(m) = x yp(m) = y Call RK4(x, y, h) m = m + 1 xp(m) = x yp(m) = y For i = 2 To n Call NSSHeun(xp(i - 2), yp(i - 2), xp(i - 1), yp(i - 1), x, y, h, iter) m = m + 1 xp(m) = x yp(m) = y iterp(m) = iter Next i Sheets("NSS Heun").Select Range("a5").Select For i = 0 To m ActiveCell.Value = xp(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = yp(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = iterp(i) ActiveCell.Offset(1, -2).Select Next i Range("a5").Select End Sub Sub RK4(x, y, h) 'Implement RK4 method Dim k1 As Single, k2 As Single, k3 As Single, k4 As Single Dim ym As Single, ye As Single, slope As Single Call Derivs(x, y, k1) ym = y + k1 * h / 2 Call Derivs(x + h / 2, ym, k2) ym = y + k2 * h / 2 Call Derivs(x + h / 2, ym, k3) ye = y + k3 * h Call Derivs(x + h, ye, k4) slope = (k1 + 2 * (k2 + k3) + k4) / 6 y = y + slope * h x = x + h End Sub Sub NSSHeun(x0, y0, x1, y1, x, y, h, iter) 'Implement Non Self-Starting Heun Dim i As Integer Dim y2 As Single Dim slope As Single, k1 As Single, k2 As Single Dim ea As Single Dim y2p As Single Static y2old As Single, y2pold As Single Call Derivs(x1, y1, k1) y2 = y0 + k1 * 2 * h y2p = y2 If iter > 0 Then y2 = y2 + 4 * (y2old - y2pold) / 5 End If x = x + h iter = 0 Do y2old = y2 Call Derivs(x, y2, k2) slope = (k1 + k2) / 2 y2 = y1 + slope * h iter = iter + 1 ea = Abs((y2 - y2old) / y2) * 100 If ea < 0.01 Then Exit Do Loop y = y2 - (y2 - y2p) / 5 y2old = y2 y2pold = y2p End Sub Sub Derivs(x, y, dydx) 'Define ODE dydx = 4 * Exp(0.8 * x) - 0.5 * y End Sub 26.12 26.13 Use Matlab to solve tspan=[0,5]'; x0=[0,0.5]'; [t,x]=ode45('dxdt',tspan,x0); plot(t,x(:,1),t,x(:,2),'--') grid title('Angle Theta and Angular Velocity Versus Time') xlabel('Time, t') ylabel('Theta (Solid) and Angular Velocity (Dashed)') axis([0 2 0 10]) zoom function dx=dxdt(t,x) dx=[x(2);(9.8/0.5)*x(1)]; A ngle Theta and A ngular V eloc ity V ers us Tim e 10 Theta (S olid) and A ngular V eloc ity (D as hed) 9 8 7 6 5 4 3 2 1 0 0 0.2 0.4 0.6 0.8 1 Tim e, t 1.2 1.4 1.6 1.8 2 A ngle Theta and A ngular V eloc ity V ers us Tim e 5 Theta (S olid) and A ngular V elo c ity (Das hed) 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 Tim e, t 26.14 Analytic solution: x = 6.004e −500t − 2 .004 e −t t=[0:.01:.02]; x=6.004*exp(-500*t)-2.004*exp(-t); plot(t,x) grid xlabel('t') ylabel('x') title('Analytic Solution:Fast Transient') gtext('6.004e^-500t-2.004 e^-t') t=[0.02:.01:5]; x=6.004*exp(-500*t)-2.004*exp(-t); plot(t,x) grid xlabel('t') ylabel('x') title('Analytic Solution: Slow Transition') gtext('6.004e^-500t-2.004 e^-t') 1.4 1.5 A naly tic S olution: F as t Trans ient 4 3 2 x 6.004e - 500t-2.004 e - t 1 0 -1 -2 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 4.5 5 t A naly tic S olution: S low Trans ition 0 -0.2 -0.4 -0.6 x -0.8 6.004e - 500t-2.004 e - t -1 -1.2 -1.4 -1.6 -1.8 -2 0 0.5 1 1.5 2 2.5 3 3.5 4 t Numerical solution: tspan=[0 5]; xo=[4]; [t,x]=ode23s('eqn',tspan,xo); plot(t,x) grid xlabel('t') ylabel('x') title('Numerical Solution: Fast Transient') axis([0 .02 -2 4]) tspan=[0 5]; xo=[4]; [t,x]=ode23s('eqn',tspan,xo); plot(t,x) grid xlabel('t') ylabel('x') title('Numerical Solution: Slow Transition') axis([0.02 5 -2 0]) Num eric al S olution: F as t Trans ient 4 3 x 2 1 0 -1 -2 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 t Num e ric al S olution: S low Trans ition 0 -0. 2 -0. 4 -0. 6 x -0. 8 -1 -1. 2 -1. 4 -1. 6 -1. 8 -2 0.5 1 1.5 2 2. 5 t 3 3.5 4 4. 5 5