Chapter 25 The Analytical Solution Can Be Derived by Separation of Variables

CHAPTER 25 25.1 The analytical solution can be derived by separation of variables ∫ dy = x 2 − 1.2 dx y ∫ ln y = x3 − 1.2 x + C 3 Substituting the initial conditions yields C = 0. Taking the exponential give the final result y= x3 −1.2 x e3 The result can be plotted as 2 1 0 0 1 2 25.2 Euler’s method with h = 0.5 x 0 0.5 1 1.5 2 y 1 0.4 0.21 0.189 0.288225 dy/dx -1.2 -0.38 -0.042 0.19845 0.80703 Euler’s method with h = 0.25 gives x 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 y 1 0.7 0.500938 0.381965 0.321089 0.305035 0.332679 0.420007 0.615572 dy/dx -1.2 -0.79625 -0.47589 -0.2435 -0.06422 0.110575 0.349312 0.782262 1.723602 The results can be plotted along with the analytical solution as 2 1 0 0 1 2 25.3 For Heun’s method, the value of the slope at x = 0 can be computed as −0.6 which can be used to compute the value of y at the end of the interval as y( 0.5) = 1 + (0 − 1.2 (1))0.5 = 0.4 The slope at the end of the interval can be computed as y' (0.5) = 0.4 (0.5) 2 − 1.2( 0.4 ) = −0.38 which can be averaged with the initial slope to predict y( 0.5) = 1 + − 0.6 − 0.38 0.5 = 0.605 2 This formula can then be iterated to yield j 0 1 2 3 4 yij εa 0.4 0.605 0.5563124 0.5678757 0.5651295 33.9 8.75 2.036 0.4859 The remaining steps can be implemented with the result xi 0.5 1.0 1.5 2.0 yi 0.5651295 0.4104059 0.5279021 2.181574 The results along with the analytical solution are displayed below: 2 1 0 0 1 2 25.4 The midpoint method with h = 0.5 x 0 0.5 1 1.5 2 y 1 0.601875 0.455591 0.534038 1.161909 dy/dx -1.2 -0.57178 -0.09112 0.56074 3.253344 ym 0.7 0.45893 0.432812 0.674223 1.975245 dy/dx-mid -0.79625 -0.29257 0.156894 1.255741 7.629383 dy/dx -1.2 -0.85121 -0.5494 -0.30293 -0.08554 0.157521 0.538615 1.344657 3.464206 ym 0.85 0.641919 0.509638 0.437323 0.417007 0.454231 0.580294 0.890046 1.670242 dy/dx-mid -1.00672 -0.68003 -0.41249 -0.18996 0.027366 0.313703 0.835986 2.061012 5.537897 with h = 0.25 gives x 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 y 1 0.74832 0.578312 0.47519 0.4277 0.434541 0.512967 0.721963 1.237216 The results can be plotted along with the analytical solution as 2 1 0 0 1 2 25.5 The 4th-order RK method with h = 0.5 gives x y k1 ym k2 0 1 -1.2 0.7 -0.79625 0.5 0.572344 -0.54373 0.436412 -0.27821 1 0.420375 -0.08407 0.399356 0.144767 1.5 0.509104 0.534559 0.642744 1.197111 2 1.29855 3.635941 2.207535 8.526606 ym 0.800938 0.50279 0.456567 0.808382 3.430202 k3 -0.91107 -0.32053 0.165505 1.505611 13.24915 ye 0.544467 0.412079 0.503128 1.26191 7.923127 k4 -0.51724 -0.08242 0.528284 3.533348 40.01179 phi -0.85531 -0.30394 0.177459 1.578892 14.53321 2 1 0 0 1 2 25.6 (a) The analytical solution can be derived by separation of variables ∫ dy y ∫ = 1 + x dx 2 y =x+ x2 +C 2 Substituting the initial conditions yields C = 2. Substituting this value and solving for y gives the final result y= ( x 2 + 2 x + 4) 2 16 The result can be plotted as 4 2 0 0 0.5 1 (b) Euler’s method with h = 0.5 x y dy/dx 0 0.5 1 1 1.5 2.418559 1 1.837117 3.110343 Euler’s method with h = 0.25 gives x 0 0.25 0.5 0.75 1 y 1 1.25 1.599386 2.073636 2.703642 dy/dx 1 1.397542 1.897002 2.520022 3.288551 The results can be plotted along with the analytical solution as 4 2 0 0 0.5 1 (c) For Heun’s method, the first step along with the associated iterations is j yij εa 0 1 2 3 1.500000 1.709279 1.740273 1.744698 12.243720 1.780954 2.536284E-01 The remaining steps can be implemented with the result xi 0.00E+00 5.00E-01 1 yi 1 1.744698 3.122586 The results along with the analytical solution are displayed below: 4 2 0 0.0 0.5 1.0 (d) The midpoint method with h = 0.5 x 0 0.5 1 y 1 1.698771 2.992924 dy/dx 1 1.955054 3.460014 ym dy/dx-mid 1.25 1.397542 2.187535 2.588305 3.857927 4.419362 dy/dx 1 1.424293 1.964934 2.645318 3.48888 ym dy/dx-mid 1.125 1.193243 1.476347 1.670694 1.961601 2.275929 2.615631 3.032421 3.479182 3.96367 with h = 0.25 gives x 0 0.25 0.5 0.75 1 y 1 1.298311 1.715984 2.284966 3.043072 The results can be plotted along with the analytical solution as 4 2 0 0 0.5 1 (e) The 4th-order RK method with h = 0.5 gives x y k1 ym k2 ym k3 ye k4 phi 0 1 1 1.25 1.397542 1.349386 1.452038 1.726019 1.970671 1.444972 0.5 1.722486 1.968653 2.214649 2.604297 2.37356 2.696114 3.070543 3.504593 2.679011 1 3.061992 3.499709 3.936919 4.464376 4.178086 4.599082 5.361533 5.788746 4.569229 4 2 0 0 0.5 1 25.7 The second-order ODE is transformed into a pair of first-order ODEs as in dy =z dx dz = x− y dx y (0)=2 z( 0)= 0 (a) The first few steps of Euler’s method are x 0 0.1 0.2 0.3 0.4 0.5 y z dy/dx dz/dx 2 0 0 -2 2 -0.2 -0.2 -1.9 1.98 -0.39 -0.39 -1.78 1.941 -0.568 -0.568 -1.641 1.8842 -0.7321 -0.7321 -1.4842 1.81099 -0.88052 -0.88052 -1.31099 (b) For Heun (without iterating the corrector) the first few steps are x y z 0 2 0 0.1 1.99 -0.195 0.2 1.96105 -0.37803 0.3 1.914442 -0.54724 0.4 1.851646 -0.70095 0.5 1.774293 -0.83761 dy/dx 0 -0.195 -0.37803 -0.54724 -0.70095 -0.83761 dz/dx -2 -1.89 -1.76105 -1.61444 -1.45165 -1.27429 yend 2 1.9705 1.923248 1.859718 1.781551 1.690532 zend -0.2 -0.384 -0.55413 -0.70868 -0.84611 -0.96504 Both results are plotted below: 4 Heun y 2 x 0 0 -2 2 z 4 Euler 25.8 The second-order ODE is transformed into a pair of first-order ODEs as in dy/dx -0.2 -0.384 -0.55413 -0.70868 -0.84611 -0.96504 dz/dx ave slope -1.9 -0.1 -1.7705 -0.2895 -1.62325 -0.46608 -1.45972 -0.62796 -1.28155 -0.77353 -1.09053 -0.90132 dy =z dt dz = −0.5z − 5 y dt y (0)=4 z (0)=0 The results for the 4th-order RK method are tabulated and plotted below: t y z 0 4.0000 0.0000 k11 k12 ymid zmid k21 k22 ymid zmid k31 k32 yend zend k41 k42 phi1 phi2 0.00 -20.00 4.00 -5.00 -5.00 -17.50 2.75 -4.38 -4.38 -11.56 1.81 -1.78 -1.78 -8.17 -3.42 -14.38 0.5 2.2891 -7.1914 -7.19 -7.85 1 -0.6167 -6.6682 -6.67 3.33 -1.04 3.95 3.95 3.23 -5.81 6.42 -2.28 -5.06 -5.06 13.95 -1.88 -3.18 -3.18 11.00 -2.21 0.49 -9.15 -9.15 4.89 4.89 8.59 -3.05 10.82 1.82 1.82 1.5 -2.1393 -1.2584 -1.26 11.33 -2.45 2.12 0.00 -6.66 -6.66 1.57 1.57 11.48 -1.75 1.61 1.61 5.75 1.16 4.93 4.93 3.67 3.67 -0.19 0.28 -1.66 -1.66 -0.55 3.16 1.23 3.49 3.49 -6.85 0.89 2.30 2.30 -5.60 1.17 -2.78 -2.78 -4.45 2.14 -5.24 -6.12 1.43 -0.14 -0.14 -7.10 1.05 -0.38 -0.38 -5.06 0.89 -1.45 -1.45 -3.75 -0.18 -5.70 -4.24 0.63 -2.52 -2.52 -1.89 0.37 -1.93 -1.93 -0.86 0.03 0.56 0.56 -0.43 -1.63 -1.70 2 -1.5614 3.3995 3.40 6.11 -0.71 2.5 0.0172 4.0139 4.01 -2.09 1.02 3 1.0852 1.3939 1.39 3.5 0.9945 -1.4562 -1.46 1.09 -0.33 7.92 -1.33 1.05 9.32 4 0.1790 -2.3048 -2.30 0.26 -0.40 -2.24 -2.24 3.11 -0.38 -1.53 -1.53 2.67 -0.59 1.51 1.51 2.17 -1.39 2.33 4.5 -0.5150 -1.1399 -1.14 3.15 -0.80 -0.35 -0.35 4.18 -0.60 -0.10 -0.10 3.07 -0.56 1.02 1.02 2.31 -0.17 3.32 2.74 -0.47 1.75 -0.30 1.01 -0.12 -0.09 -0.09 0.65 1.49 5 -0.6001 0.5213 0.52 1.21 1.21 0.96 0.96 8 y 4 t 0 0 2 4 -4 z -8 25.9 (a) The Heun method without iteration can be implemented as in the following table: t 0 0.1 0.2 0.3 0.4 0.5 • • • 2.9 3 y 1 1.000498 1.002971 1.00933 1.021391 1.040874 k1 0 0.009972 0.039587 0.088147 0.15489 0.239244 yend 1 1.001496 1.00693 1.018145 1.03688 1.064798 k2 0.009967 0.039489 0.087592 0.153062 0.234765 0.331852 4.527257 4.544722 0.259141 0.090507 4.553171 4.553773 0.09016 0.007858 phi 0.004983 0.02473 0.063589 0.120604 0.194828 0.285548 • • • 0.17465 0.049183 (b) The Ralston 2nd order RK method can be implemented as in the following table: t y k1 yint k2 phi 0.79 0 0.1 0.2 0.3 0.4 0.5 • • • 2.9 3 1 1.000374 1.00273 1.008993 1.021013 1.040573 0 0.00997 0.039577 0.088118 0.154833 0.239175 1 1.001122 1.005698 1.015602 1.032626 1.058511 0.005614 0.030348 0.074158 0.136249 0.215982 0.313061 4.779856 4.797775 0.2736 0.095547 4.800376 4.804941 0.131997 0.021276 0.003743 0.023555 0.062631 0.120205 0.195599 0.288432 • • • 0.179198 0.046033 Both methods are displayed on the following plot along with the exact solution. The Ralston method performs much better for this case. Ralston 4 Heun 2 0 0 1 2 3 25.10 The solution results are as in the following table and plot: t 0 0.5 1 1.5 2 2.5 3 y 1 0.640625 0.626628 1.019213 1.855358 3.157404 4.938848 k1 -1 -0.39063 0.373372 1.230787 2.144642 3.092596 4.061152 k2 -0.6875 0.019531 0.842529 1.735591 2.670982 3.631947 4.608364 k3 -0.5625 0.144531 0.967529 1.860591 2.795982 3.756947 4.733364 phi -0.71875 -0.02799 0.78517 1.67229 2.604092 3.562889 4.537995 4 2 0 0 1 2 25.11 (a) Euler x 0 y 2.0000 z 4.0000 dy/dx 1.00 dz/dx -16.00 3 0.2 0.4 0.6 0.8 1 2.2000 2.1387 1.9536 1.7209 1.4819 0.8000 0.6592 0.5663 0.5036 0.4600 -0.31 -0.93 -1.16 -1.20 -1.12 -0.70 -0.46 -0.31 -0.22 -0.16 (b) 4th-order RK x y 0 0.2 0.4 0.6 0.8 1 z 2.000 2.083 2.003 1.840 1.641 1.433 4.000 2.278 1.571 1.168 0.923 0.768 k11 k12 1.000 -16.000 -0.071 -5.406 -0.655 -2.472 -0.937 -1.256 -1.035 -0.699 -1.027 -0.423 k21 0.324 -0.447 -0.843 -1.010 -1.042 -0.997 k22 -6.048 -3.134 -1.698 -0.950 -0.559 -0.351 k31 k32 0.459 -11.714 -0.372 -3.934 -0.806 -1.884 -0.996 -1.002 -1.040 -0.577 -1.003 -0.358 k41 -0.090 -0.665 -0.941 -1.036 -1.026 -0.960 k42 -0.123 -1.686 -2.438 -2.207 -1.667 -1.143 phi1 0.413 -0.396 -0.816 -0.997 -1.038 -0.998 phi2 -8.608 -3.538 -2.012 -1.228 -0.773 -0.497 Both methods are plotted on the same graph below. Notice how Euler’s method (particularly for z) is very inaccurate for this step size. The 4th-order RK is much closer to the exact solution. 4 4th-order RK Euler y 2 z 0 0 0.5 ( x − 2 )2 − 25.12 dy = 10e 2 ( 0.075)2 − 0.6 y dx 4th-order RK method: One step (h = 0.5): Two steps (h = 0.25): y1 = 0.3704188 y2 = 0.3704096 ∆present = −9.119×10−6 correction = ∆ = −6.08 × 10 −7 15 y2 = 0.370409 dy = −0.3 dx y scale = 0.5 + 0.5( −0.3) = 0.65 1 ∆new = 0.001(0.65) = 0.00065 Since ∆present < ∆new, therefore, increase step. hnew 0.00065 = 0.5 9.119 × 10 −6 0.2 = 11737 . 25.13 We will look at the first step only ∆present = y2 − y1 = −0.24335 dy dx = 4 e 0 − 0.5( 2) = 3 yscale = 2 + (2(3)) = 8 ∆new = 0.001(8) = 0.008 Because ∆present > ∆new, decrease step. 25.14 The calculation of the k’s can be summarized in the following table: k1 k2 k3 k4 k5 k6 x 0 0.25 0.375 0.923077 1 0.5 y 2 2.75 3.268609 5.636774 5.878223 3.805418 f(x,y) 3 3.510611 3.765131 5.552467 5.963052 4.06459 k 3 3.510611 3.765131 5.552467 5.963052 4.06459 These can then be used to compute the 4th-order prediction 1408 2197 1  25  y1 = 2 +  3+ 3.765131 + 5.552467 − 5.963052  1 = 6.193807  216  2565 4104 5 along with a fifth-order formula:  16 6656 28,561 9 2  y1 = 2 +  3+ 3.765131 + 5.552467 − 5.963052 + 4.06459  1 = 6.194339  135  12 ,825 56,430 50 55 The error estimate is obtained by subtracting these two equations to give E a = 6.194339 − 6.193807 = 0.000532 25.15 Option Explicit Sub EulerTest() Dim i As Integer, m As Integer Dim xi As Single, yi As Single, xf As Single, dx As Single, xout As Single Dim xp(200) As Single, yp(200) As Single 'Assign values yi = 1 xi = 0 xf = 4 dx = 0.5 xout = 0.5 '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 Euler(x, y, h, ynew) y = ynew If (x >= xend) Then Exit Do Loop End Sub 'Trim step if increment exceeds end Sub Euler(x, y, h, ynew) Dim dydx As Single 'Implement Euler's method Call Derivs(x, y, dydx) ynew = y + dydx * h x = x + h End Sub Sub Derivs(x, y, dydx) 'Define ODE dydx = -2 * x ^ 3 + 12 * x ^ 2 - 20 * x + 8.5 End Sub 25.16 Example 25.1: Example 25.4 (nonlinear model). Change time steps and initial conditions to 'Assign values yi = 0 xi = 0 xf = 15 dx = 0.5 xout = 1 Change Derivs Sub to Sub Derivs(t, v, dvdt) 'Define ODE dvdt = 9.8 - 12.5 / 68.1 * (v + 8.3 * (v / 46) ^ 2.2) End Sub 25.17 Option Explicit Sub Dim Dim Dim RK4Test() 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 = 1 xi = 0 xf = 4 dx = 0.5 xout = 0.5 '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 RK4(x, y, h, ynew) y = ynew If (x >= xend) Then Exit Do Loop End Sub 'Trim step if increment exceeds end Sub RK4(x, y, h, ynew) '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 ynew = y + slope * h x = x + h End Sub Sub Derivs(x, y, dydx) 'Define ODE dydx = -2 * x ^ 3 + 12 * x ^ 2 - 20 * x + 8.5 End Sub 25.18 Example 25.1: Example 25.5 Change time steps and initial conditions to 'Assign values yi = 2 xi = 0 xf = 4 dx = 1 xout = 1 Change Derivs Sub to Sub Derivs(x, y, dydx) 'Define ODE dydx = 4 * Exp(0.8 * x) - 0.5 * y End Sub 25.19 Option Explicit Sub Dim Dim Dim RK4SysTest() 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 = 2 yi(1) = 4 yi(2) = 6 dx = 0.5 xout = 0.5 '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 RK4Sys(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 RK4Sys(x, y, h, n, ynew) Dim j As Integer Dim dydx(10) As Single Dim ym(10), ye(10) Dim k1(10), k2(10), k3(10), k4(10) Dim slope(10) 'Implement RK4 method for systems of ODEs Call Derivs(x, y, k1()) For j = 1 To n ym(j) = y(j) + k1(j) * h / 2 Next j Call Derivs(x + h / 2, ym, k2()) For j = 1 To n ym(j) = y(j) + k2(j) * h / 2 Next j Call Derivs(x + h / 2, ym, k3()) For j = 1 To n ye(j) = y(j) + k3(j) * h Next j Call Derivs(x + h, ye, k4()) For j = 1 To n slope(j) = (k1(j) + 2 * (k2(j) + k3(j)) + k4(j)) / 6 Next j For j = 1 To n ynew(j) = y(j) + slope(j) * h Next j x = x + h End Sub Sub Derivs(x, y, dydx) 'Define ODE dydx(1) = -0.5 * y(1) dydx(2) = 4 - 0.3 * y(2) - 0.1 * y(1) End Sub Application to Example 25.10: 25.20 Main Program: %Damped spring mass system %mass: m=10 kg %damping: c=5,40,200 N/(m/s) %spring: k=40 N/m % MATLAB 5 version %Independent Variable t, tspan=[tstart tstop] %initial conditions [x(1)=velocity, x(2)=displacement]; tspan=[0 15]; ic=[0 1]; global cm km m=10; c(1)=5; c(2)=40; c(3)=200; k=40; km=k/m; end for n=1:3 cm=c(n)/m [t,x]=ode45('kc',tspan,ic); plot(t,x(:,2)); grid; xlabel('time - sec.'); ylabel('displacement - m'); title('m(d2x/dt2)+c(dx/dt)+kx=0; m=10 kg, k= 40 N/m') hold on gtext('c=5');gtext('cc=40');gtext('c=200 N/(m/s)') Function ‘kc’: %Damped spring mass system - m d2x/dt2 + c dx/dt + k x =0 %mass: m=10 kg %damping: c=5,40,200 N/(m/s) %spring: k=40 N/m %x(1)=velocity, x(2)=dispacement function dx=kc(t,x); global cm km dx=[-cm*x(1)-km*x(2); x(1)]; m (d2x /dt2)+ c (dx /dt)+ k x = 0; m = 10 k g, k = 40 N/m 1 0.8 c = 200 N/(m /s ) 0.6 c c = 40 dis plac em ent - m 0.4 0.2 0 -0.2 -0.4 c=5 -0.6 -0.8 0 5 10 15 t im e - s ec . 25.21 The Matlab program on the following pages performs the Euler Method and the plots shows the depth of the water vs. time. From the plot, we approximate that it takes about 58 minutes to drain the cylindrical tank. %euler.m dt=0.5; max=60; n=max/dt+1; t=zeros(1,n); y=zeros(1,n); t(1)=0; y(1)=9; for i=1:n y(i+1)=y(i)+dydt(t(i),y(i))*dt; t(i+1)=t(i)+dt; end plot(t,y) grid xlabel('Time-minutes') ylabel('Depth of the Water-ft') title('How Long Does it Take to Drain a 9-ft high Cylindrical Tank?') zoom function dy=dydt(t,y); dy=-0.1*sqrt(y); 25.22 x = x(1) dx v= = x (2 ) dt dv d  dx  d 2 x =  = dt dt  dt  dt 2 dx(1) = x (2) dt dx( 2) = −5 x(1) x (2) − ( x (1) + 7 ) sin(t ) dt x (1)(t = 0) = 6 x (2)(t = 0) = 1.5 tspan=[0,15]'; x0=[6,1.5]'; [t,x]=ode45('dxdt',tspan,x0); plot(t,x(:,1),t,x(:,2),'--') grid title('Displacement and Velocity Versus Time') xlabel('Time, t') ylabel('Displacement and Velocity') gtext('displacement') gtext('velocity') function dx=dxdt(t,x) dx=[x(2);-5*x(1)*x(2)+(x(1)+7)*sin(1*t)]; D is plac em ent and V eloc ity V ers us Tim e 7 6 dis plac em ent Dis plac em ent and V eloc ity 5 4 3 2 1 veloc ity 0 -1 0 5 10 Tim e, t 15