Chapter 27 The Solution Can Be Assumed to Be T = elx

CHAPTER 27 27.1 The solution can be assumed to be T = eλx. This, along with the second derivative T” = λ2eλx, can be substituted into the differential equation to give λ2 e λx − 0.1e λx = 0 which can be used to solve for λ2 − 0.1 = 0 λ = ± 0.1 Therefore, the general solution is T = Ae 0.1x + Be − 0.1x The constants can be evaluated by substituting each of the boundary conditions to generate two equations with two unknowns, 200 = A + B 100 = 23.62434 A + 0.042329 B which can be solved for A = 3.881524 and B = 196.1185. The final solution is, therefore, T = 3.881524 e 0.1x + 196.1185e − 0.1x which can be used to generate the values below: x 0 1 2 3 4 5 6 7 8 9 10 T 200 148.2747 111.5008 85.97028 69.10864 59.21565 55.29373 56.94741 64.34346 78.22764 100 200 100 0 0 2 4 6 8 10 27.2 Reexpress the second-order equation as a pair of ODEs: dT =z dx dz = 0.1T dx The solution was then generated on the Excel spreadsheet using the Heun method (without iteration) with a step-size of 0.01. An initial condition of z = −55 was chosen for the first shot. The first few calculation results are shown below. x 0 0.1 0.2 0.3 0.4 0.5 T z 200.000 194.600 189.395 184.378 179.547 174.894 -55.000 -53.028 -51.108 -49.240 -47.420 -45.649 k11 k12 Tend zend k21 k22 phi1 phi2 -55.000 20.000 194.500 -53.000 -53.000 19.450 -54.000 19.725 -53.028 19.460 189.297 -51.082 -51.082 18.930 -52.055 19.195 -51.108 18.939 184.284 -49.214 -49.214 18.428 -50.161 18.684 -49.240 18.438 179.454 -47.396 -47.396 17.945 -48.318 18.192 -47.420 17.955 174.805 -45.625 -45.625 17.480 -46.523 17.718 -45.649 17.489 170.330 -43.900 -43.900 17.033 -44.774 17.261 The resulting value at x = 10 was T(10) = 315.759. A second shot using an initial condition of z (0) = −70 was attempted with the result at x = 10 of T(10) = −243.249. These values can then be used to derive the correct initial condition, z( 0) = −55 + − 70 + 55 (100 − 315.759 ) = −60.79 − 243.249 − 315.759 The resulting fit, along with the two “shots” are displayed below: 400 300 200 100 0 -100 0 2 4 6 8 10 -200 -300 27.3 A centered finite difference can be substituted for the second derivative to give, Ti −1 − 2Ti + Ti +1 h2 − 0.1Ti = 0 or for h = 1, − Ti −1 + 2.1Ti − Ti +1 = 0 The first node would be 2.1T1 − T2 = 200 and the last node would be − T9 + 2.1T10 = 100 The tridiagonal system can be solved with the Thomas algorithm or Gauss-Seidel for (the analytical solution is also included) x 0 1 2 3 4 5 6 7 8 9 10 T 200 148.4838 111.816 86.32978 69.47655 59.57097 55.62249 57.23625 64.57365 78.3684 100 Analytical 200 148.2747 111.5008 85.97028 69.10864 59.21565 55.29373 56.94741 64.34346 78.22764 100 27.4 The second-order ODE can be expressed as the following pair of first-order ODEs, dy =z dx dz 2 z + y − x = dx 8 These can be solved for two guesses for the initial condition of z. For our cases we used z(0) −1 −0.5 y(20) −6523.000507 7935.937904 Clearly, the solution is quite sensitive to the initial conditions. These values can then be used to derive the correct initial condition, z( 0) = −1 + − 0.5 + 1 (8 + 6523.000507 ) = −0.774154 7935.937904 + 6523.000507 The resulting fit is displayed below: 12 8 4 0 0 5 10 15 20 27.5 Centered finite differences can be substituted for the second and first derivatives to give, 8 yi +1 − 2 y i + yi −1 ∆x 2 −2 y i +1 − y i −1 − yi + x i = 0 ∆x or substituting ∆x = 2 and collecting terms yields 2.5 y i +1 − 5 yi + 1.5 yi −1 + x i = 0 This equation can be written for each node and solved with either the Gauss-Seidel method or a tridiagonal solver to give x 0 2 4 6 8 10 12 14 16 18 20 T 5 4.287065 4.623551 5.600062 6.960955 8.536414 10.18645 11.72749 12.78088 12.39044 8 12 8 4 0 0 5 10 15 20 27.6 The second-order ODE can be expressed as the following pair of first-order ODEs, dT =z dx dz 7 4 = 1.2 × 10 (T + 273) − 5(150 − T ) dx The solution was then generated on the Excel spreadsheet using the Heun method (without iteration) with a step-size of 0.01. The Excel Solver was used to adjust the initial condition of z until the value of T(0.5) = 100. Part of the resulting spreadsheet is shown below along with a graph of the final solution. x T z k11 k12 Tend zend k21 k22 φ1 φ2 0 200.000 -927.673 -927.673 6256.560 190.723 -865.107 -865.107 5752.643 -896.390 6004.601 0.01 191.036 -867.627 -867.627 5769.196 182.360 -809.935 -809.935 5321.210 -838.781 5545.203 0.02 0.03 0.04 0.05 182.648 174.793 167.433 160.532 -812.175 -760.816 -713.115 -668.694 -812.175 -760.816 -713.115 -668.694 5335.738 4948.905 4602.594 4291.667 174.527 167.185 160.301 153.845 -758.817 -711.327 -667.089 -625.778 -758.817 -711.327 -667.089 -625.778 4936.083 4591.217 4281.522 4002.685 -785.496 -736.071 -690.102 -647.236 5135.910 4770.061 4442.058 4147.176 200 100 0 0 0.1 0.2 0.3 0.4 0.5 27.7 The second-order ODE can be linearized as in 2 d T 7 4 7 3 − 1.2 × 10 ( Tb + 273) − 4.8 × 10 ( Tb + 273) ( T − Tb ) + 5(150 − T ) = 0 2 dx Substituting Tb = 150 and collecting terms gives d 2T dx 2 − 41.32974 T + 2357 .591 = 0 Substituting a centered-difference approximation of the second derivative gives − Ti −1 + (2 + 41.32974 ∆x 2 ) Ti − Ti +1 = 2357.591∆x 2 We used the Gauss-Seidel method to solve these equations. The results for a few selected points are: x T 0 200 0.1 134.2765 0.2 101.5758 0.3 0.4 87.91595 87.45616 0.5 100 A graph of the entire solution along with the nonlinear result from Prob. 27.7 is shown below: 200 Linear 100 Nonlinear 0 0 0.1 0.2 0.3 0.4 0.5 27.8 For three springs  2k1  − ω 2  A1   m1  k − 1 A1 m1 − k1 m1 = 0 A2  2k  +  1 − ω 2  A2  m1  k − 1 A2 m1 − k1 A3 m1 = 0  2k  +  1 − ω 2  A3  m1  = 0 Substituting m = 40 kg and k = 240 gives ( 12 − ω ) A 2 − 6 A1 1 ( − 6 A2 ) + 12 − ω 2 A2 − 6 A2 ( − 6 A3 + 12 − ω 2 )A 3 = 0 = 0 = 0 The determinant is − ω 6 + 36ω 4 − 360ω 2 + 864 = 0 which can be solved for ω2 = 20.4853, 12, and 3.5147 s−2. Therefore the frequencies are ω = 4.526, 3.464, and 1.875 s−1. Substituting these values into the original equations yields for ω2 = 20.4853, A1 = −0.707A2 = A3 for ω2 = 12 A1 = −A3, and = A2 = 0 for ω2 = 3.5147 A1 = 0.707A2 = A3 Plots: 0 0 0 4 4 4 27.9 For 5 interior points (h = 3/6 = 0.5), the result is Eq. (27.19) with 2 − 0.25p2 on the diagonal. Dividing by 0.25 gives, 8 − p 2  −4      −4 8 − p2 −4 −4 8 − p2 −4 −4 8 − p2 −4    =0 −4  2 8− p  The determinant can be expanded (e.g., with Fadeev-Leverrier or the MATLAB poly function) to give 0 = − p 10 + 40 p 8 − 576 p 6 + 3,584 p 4 − 8960 p 2 + 6,144 The roots of this polynomial can be determined as (e.g., with Bairstow’s methods or the MATLAB roots function) p2 = 1.072, 4, 8, 12, 14.94. The square root of these roots yields p = 1.035, 2, 2.828, 3.464, and 3.864. 27.10 Minors: (2 − λ ) 3 −4 λ 4 8 4 8 3 − λ = −λ 3 + 10λ 2 + 101λ + 18 7 − λ − 2 10 7 − λ + 10 10 4 27.11 Although the following computation can be implemented on a pocket calculator, a spreadsheet or with a program, we’ve used MATLAB. >> a=[2 2 10;8 3 4;10 4 5] a = 2 2 10 8 3 4 10 4 5 >> x=[1 1 1]' x = 1 1 1 First iteration: >> x=a*x x = 14 15 19 >> e=max(x) e = 19 >> x=x/e x = 0.7368 0.7895 1.0000 Second iteration: >> x=a*x x = 13.0526 12.2632 15.5263 >> e=max(x) e = 15.5263 >> x=x/e x = 0.8407 0.7898 1.0000 Third iteration: >> x=a*x x = 13.2610 13.0949 16.5661 >> e=max(x) e = 16.5661 >> x=x/e x = 0.8005 0.7905 1.0000 Fourth iteration: >> x=a*x x = 13.1819 12.7753 16.1668 >> e=max(x) e = 16.1668 >> x=x/e x = 0.8154 0.7902 1.0000 Thus, after four iterations, the result is converging on a highest eigenvalue of 16.2741 with a corresponding eigenvector of [0.811 0.790 1]. 27.12 As in Example 27.10, the computation can be laid out as 2 8 10 2 3 4 10 4 5 First iteration: -0.05556 1.666667 0 -5 0.111111 0.666667 -1.22222 4 -0.55556 1 1 1 Second iteration: -0.05556 1.666667 0 -5 -1.22222 4 -0.38889 1 [A] = = = 0.388889 -1 0.222222 1.959877 -5.88889 eigenvalue eigenvector -0.3888889 -1 1 -0.2222222 -5.88889 -0.3328092 1 0.111111 0.666667 -0.55556 -0.22222 0.746914 -0.1268344 Third iteration: -0.05556 1.666667 0 -5 0.111111 0.666667 -1.22222 4 -0.55556 -0.33281 1 -0.12683 = 1.840176 -5.50734 0.700151 -5.50734 -0.3341317 1 -0.1271307 Fourth iteration: -0.05556 1.666667 0 -5 0.111111 0.666667 -1.22222 4 -0.55556 -0.33413 1 -0.12713 = 1.840611 -5.50852 0.700169 -5.50852 -0.3341389 1 -0.1271065 Thus, after four iterations, the estimate of the lowest eigenvalue is 1/(−5.5085) = −0.1815 with an eigenvector of [−0.3341 1 -0.1271]. 27.13 Here is VBA Code to implement the shooting method: Public hp As Single, Ta As Single Option Explicit Sub Shoot() Dim Dim Dim Dim Dim Dim Dim n As Integer, m As Integer, i As Integer, j As Integer x0 As Single, xf As Single x As Single, y(2) As Single, h As Single, dx As Single, xend As Single xp(200) As Single, yp(2, 200) As Single, xout As Single z01 As Single, z02 As Single, T01 As Single, T02 As Single T0 As Single, Tf As Single Tf1 As Single, Tf2 As Single 'set parameters n = 2 hp = 0.01 Ta = 20 x0 = 0 T0 = 40 xf = 10 Tf = 200 dx = 2 xend = xf xout = 2 'first shot x = x0 y(1) = T0 y(2) = 10 Call RKsystems(x, y, n, dx, xf, z01 = yp(2, 0) Tf1 = yp(1, m) 'second shot x = x0 y(1) = T0 y(2) = 20 Call RKsystems(x, y, n, dx, xf, z02 = yp(2, 0) Tf2 = yp(1, m) 'last shot x = x0 y(1) = T0 'linear interpolation y(2) = z01 + (z02 - z01) / (Tf2 Call RKsystems(x, y, n, dx, xf, 'output results Range("a4:C1004").ClearContents Range("A4").Select For j = 0 To m ActiveCell.Value = xp(j) xout, xp(), yp(), m) xout, xp(), yp(), m) - Tf1) * (Tf - Tf1) xout, xp(), yp(), m) For i = 1 To n ActiveCell.Offset(0, 1).Select ActiveCell.Value = yp(i, j) Next i ActiveCell.Offset(1, -n).Select Next j Range("A4").Select End Sub Sub RKsystems(x, y, n, dx, xf, xout, xp, yp, m) Dim i As Integer Dim xend As Single, h As Single m = 0 For i = 1 To n yp(i, m) = y(i) Next i Do xend = x + xout If xend > xf Then xend = xf h = dx Do If xend - x < h Then h = xend - x Call RK4(x, y, n, h) If x >= xend Then Exit Do Loop m = m + 1 xp(m) = x For i = 1 To n yp(i, m) = y(i) Next i If x >= xf Then Exit Do Loop End Sub Sub RK4(x, y, n, h) Dim i Dim ynew, dydx(10), ym(10), ye(10) Dim k1(10), k2(10), k3(10), k4(10) Dim slope(10) Call Derivs(x, y, k1) For i = 1 To n ym(i) = y(i) + k1(i) * h / 2 Next i Call Derivs(x + h / 2, ym, k2) For i = 1 To n ym(i) = y(i) + k2(i) * h / 2 Next i Call Derivs(x + h / 2, ym, k3) For i = 1 To n ye(i) = y(i) + k3(i) * h Next i Call Derivs(x + h, ye, k4) For i = 1 To n slope(i) = (k1(i) + 2 * (k2(i) + k3(i)) + k4(i)) / 6 Next i For i = 1 To n y(i) = y(i) + slope(i) * h Next i x = x + h End Sub Sub Derivs(x, y, dydx) dydx(1) = y(2) dydx(2) = hp * (y(1) - Ta) End Sub 27.14 27.15 A general formulation that describes Example 27.3 as well as Probs. 27.3 and 27.5 is d2 y dy a 2 + b + cy + f ( x) = 0 dx dx Finite difference approximations can be substituted for the derivatives: a yi +1 − 2 y i + y i −1 y − y i −1 + b i +1 + cyi + f ( xi ) = 0 2 ∆x 2∆x Collecting terms ( ) − ( a − 0.5b∆x ) y i −1 + 2a + c∆x 2 y i − ( a + 0.5b∆x ) y i +1 = f ( xi )∆x 2 The following VBA code implants this equation as applied to Example 27.3. Public hp As Single, dx As Single Option Explicit Sub FDBoundaryValue() Dim ns As Integer, i As Integer Dim a As Single, b As Single, c As Single Dim e(100) As Single, f(100) As Single, g(100) As Single Dim r(100) As Single, y(100) As Single Dim Lx As Single, xx As Single, x(100) As Single Lx = 10 dx = 2 ns = Lx / dx xx = 0 hp = 0.01 a = 1 b = 0 c = hp y(0) = 40 y(ns) = 200 For i = 0 To ns x(i) = xx xx = xx + dx Next i f(1) = 2 * a / dx ^ 2 + c g(1) = -(a / dx ^ 2 - b / (2 * dx)) r(1) = ff(x(1)) + (a / dx ^ 2 + b / (2 * dx)) * y(0) For i = 2 To ns - 2 e(i) = -(a / dx ^ 2 + b / (2 * dx)) f(i) = 2 * a / dx ^ 2 + c g(i) = -(a / dx ^ 2 - b / (2 * dx)) r(i) = ff(x(i)) Next i e(ns - 1) = -(a / dx ^ 2 + b / (2 * dx)) f(ns - 1) = 2 * a / dx ^ 2 + c r(ns - 1) = ff(x(ns - 1)) + (a / dx ^ 2 - b / (2 * dx)) * y(ns) Sheets("Sheet2").Select Range("a5:d105").ClearContents Range("a5").Select For i = 1 To ns - 1 ActiveCell.Value = e(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = f(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = g(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = r(i) ActiveCell.Offset(1, -3).Select Next i Range("a5").Select Call Tridiag(e, f, g, r, ns - 1, y) Sheets("Sheet1").Select Range("a5:b105").ClearContents Range("a5").Select For i = 0 To ns ActiveCell.Value = x(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = y(i) ActiveCell.Offset(1, -1).Select Next i Range("a5").Select End Sub Sub Tridiag(e, f, g, r, n, x) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub Function ff(x) ff = hp * 20 End Function 27.16 27.17 Option Explicit Sub Dim Dim Dim Dim Dim Power() n As Integer, i As Integer, iter As Integer aa As Single, bb As Single a(10, 10) As Single, c(10) As Single lam As Single, lamold As Single, v(10) As Single es As Single, ea As Single es = 0.001 n = 3 aa = 2 / 0.5625 bb = -1 / 0.5625 a(1, 1) = aa a(1, 2) = bb For i = 2 To n - 1 a(i, i - 1) = bb a(i, i) = aa a(i, i + 1) = bb Next i a(i, i - 1) = bb a(i, i) = aa lam = 1 For i = 1 To n v(i) = lam Next i Sheets("sheet1").Select Range("a3:b1000").ClearContents Range("a3").Select Do iter = iter + 1 Call Mmult(a(), (v()), v(), n, n, 1) lam = Abs(v(1)) For i = 2 To n If Abs(v(i)) > lam Then lam = Abs(v(i)) Next i ActiveCell.Value = "iteration: " ActiveCell.Offset(0, 1).Select ActiveCell.Value = iter ActiveCell.Offset(1, -1).Select ActiveCell.Value = "eigenvalue: " ActiveCell.Offset(0, 1).Select ActiveCell.Value = lam ActiveCell.Offset(1, -1).Select For i = 1 To n v(i) = v(i) / lam Next i ActiveCell.Value = "eigenvector:" ActiveCell.Offset(0, 1).Select For i = 1 To n ActiveCell.Value = v(i) ActiveCell.Offset(1, 0).Select Next i ActiveCell.Offset(1, -1).Select ea = Abs((lam - lamold) / lam) * 100 lamold = lam If ea <= es Then Exit Do Loop End Sub Sub Mmult(a, b, c, m, n, l) Dim i As Integer, j As Integer, k As Integer Dim sum As Single For i = 1 To n sum = 0 For k = 1 To m sum = sum + a(i, k) * b(k) Next k c(i) = sum Next i End Sub • • • 27.18 Option Explicit Sub Dim Dim Dim Dim Dim Dim Power() n As Integer, i As Integer, iter As Integer, j As Integer aa As Single, bb As Single a(10, 10) As Single, c(10) As Single lam As Single, lamold As Single, v(10) As Single es As Single, ea As Single x(10) As Single, ai(10, 10) As Single es = 0.0000011 n = 3 aa = 2 / 0.5625 bb = -1 / 0.5625 a(1, 1) = aa a(1, 2) = bb For i = 2 To n - 1 a(i, i - 1) = bb a(i, i) = aa a(i, i + 1) = bb Next i a(i, i - 1) = bb a(i, i) = aa Call LUDminv(a(), n, x()) lam = 1 For i = 1 To n v(i) = lam Next i Sheets("sheet1").Select Range("a3:j1000").ClearContents Range("a3").Select ActiveCell.Value = "Matrix inverse:" ActiveCell.Offset(1, 0).Select For i = 1 To n For j = 1 To n ActiveCell.Value = a(i, j) ActiveCell.Offset(0, 1).Select Next j ActiveCell.Offset(1, -n).Select Next i ActiveCell.Offset(1, 0).Select Do iter = iter + 1 Call Mmult(a(), (v()), v(), n, n, 1) lam = Abs(v(1)) For i = 2 To n If Abs(v(i)) > lam Then lam = Abs(v(i)) Next i ActiveCell.Value = "iteration: " ActiveCell.Offset(0, 1).Select ActiveCell.Value = iter ActiveCell.Offset(1, -1).Select ActiveCell.Value = "eigenvalue: " ActiveCell.Offset(0, 1).Select ActiveCell.Value = lam ActiveCell.Offset(1, -1).Select For i = 1 To n v(i) = v(i) / lam Next i ActiveCell.Value = "eigenvector:" ActiveCell.Offset(0, 1).Select For i = 1 To n ActiveCell.Value = v(i) ActiveCell.Offset(1, 0).Select Next i ActiveCell.Offset(1, -1).Select ea = Abs((lam - lamold) / lam) * 100 lamold = lam If ea <= es Then Exit Do Loop End Sub Sub Mmult(a, b, c, m, n, l) Dim i As Integer, j As Integer, k As Integer Dim sum As Single For i = 1 To n sum = 0 For k = 1 To m sum = sum + a(i, k) * b(k) Next k c(i) = sum Next i End Sub Sub Dim Dim Dim LUDminv(a, n, x) i As Integer, j As Integer, er As Integer o(3) As Single, s(3) As Single, b(3) As Single ai(10, 10) As Single, tol As Single tol = 0.00001 Call Decompose(a, n, tol, o(), s(), er) If er = 0 Then For i = 1 To n For j = 1 To n If i = j Then b(j) = 1 Else b(j) = 0 End If Next j Call Substitute(a, o, n, b, x) For j = 1 To n ai(j, i) = x(j) Next j Next i End If For i = 1 To n For j = 1 To n a(i, j) = ai(i, j) Next j Next i End Sub Sub Decompose(a, n, tol, o, s, er) Dim i As Integer, j As Integer, k As Integer Dim factor As Single For i = 1 To n o(i) = i s(i) = Abs(a(i, 1)) For j = 2 To n If Abs(a(i, j)) > s(i) Then s(i) = Abs(a(i, j)) Next j Next i For k = 1 To n - 1 Call Pivot(a, o, s, n, k) If Abs(a(o(k), k) / s(o(k))) < tol Then er = -1 Exit For End If For i = k + 1 To n factor = a(o(i), k) / a(o(k), k) a(o(i), k) = factor For j = k + 1 To n a(o(i), j) = a(o(i), j) - factor * a(o(k), j) Next j Next i Next k If (Abs(a(o(k), k) / s(o(k))) < tol) Then er = -1 End Sub Sub Pivot(a, o, s, n, k) Dim ii As Integer, p As Integer Dim big As Single, dummy As Single p = k big = Abs(a(o(k), k) / s(o(k))) For ii = k + 1 To n dummy = Abs(a(o(ii), k) / s(o(ii))) If dummy > big Then big = dummy p = ii End If Next ii dummy = o(p) o(p) = o(k) o(k) = dummy End Sub Sub Substitute(a, o, n, b, x) Dim k As Integer, i As Integer, j As Integer Dim sum As Single, factor As Single For k = 1 To n - 1 For i = k + 1 To n factor = a(o(i), k) b(o(i)) = b(o(i)) - factor * b(o(k)) Next i Next k x(n) = b(o(n)) / a(o(n), n) For i = n - 1 To 1 Step -1 sum = 0 For j = i + 1 To n sum = sum + a(o(i), j) * x(j) Next j x(i) = (b(o(i)) - sum) / a(o(i), i) Next i End Sub • • • 27.19 This problem can be solved by recognizing that the solution corresponds to driving the differential equation to zero. To do this, a finite difference approximation can be substituted for the second derivative to give R= Ti −1 − 2 Ti + Ti +1 ( ∆x ) 2 − 1.2 × 10 −7 ( Ti + 273) 4 + 5(150 − Ti ) where R = the residual, which is equal to zero when the equation is satisfied. Next, a spreadsheet can be set up as below. Guesses for T can be entered in cells B11:B14. Then, the residual equation can be written in cells C11:C14 and referenced to the temperatures in column B. The square of the R’s can then be entered in column D and summed (D17). The Solver can then be invoked to drive cell D17 to zero by varying B11:B14. The result is as shown in the spreadsheet. A plot is also displayed below. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 A E sigma k Ta T0 Tn dx B C D 1 1.20E-07 5 150 200 100 0.1 x T R 0 200 0.1 133.015 4.32E-05 0.2 97.79076 0.000185 0.3 82.63883 -0.00076 0.4 83.31515 0.001114 0.5 100 SSR R^2 1.87E-09 3.42E-08 5.8E-07 1.24E-06 =sum(D11:D14) 1.86E-06 =(B10-2*B11+B12)/$B$7^2-$B$2*(B11+273)^4+$B$3*($B$4-B11) 200 150 100 50 0 0 0.1 0.2 0.3 0.4 0.5 27.20 First, an m-file containing the system of ODEs can be created and saved (in this case as odesys.m), function dy = predprey(t,y) dy=[0.3*y(1)-1.5*y(1)*y(2);-0.1*y(2)+0.036*y(1)*y(2)]; Then, the following MATLAB session is used to generate the solution: >> [t,y]=ode45('odesys',[0 100],[1;.05]); A plot of the solution along with the state-space plot are generated with >> plot(t,y) >> plot(y(:,1),y(:,2)) These plots are displayed below 14 12 0.5 10 8 6 4 2 0 0 0 20 40 60 80 100 0 10 27.21 First, the 2nd-order ODE can be reexpressed as the following system of 1st-order ODE’s dx =z dt dz = −8.333333z − 1166.667 x dt Next, we create an m-file to hold the ODEs: function dx=spring(t,y) dx=[y(2);-8.333333*y(2)-1166.667*y(1)] Then we enter the following commands into MATLAB [t,y]=ode45('spring',[0 .4],[0.3;0]) plot(t,y(:,1)); The following plot results: 0.2 0 -0.2 -0.4 0 0.1 0.2 0.3 0.4 (b) The eigenvalues and eigenvectors can be determined with the following commands: >> a=[0 -1;8.333333 1166.667]; >> format short e >> [v,d]=eig(a) v = -9.9997e-001 8.5715e-004 7.1427e-003 -1.0000e+000 d = 7.1429e-003 0 0 1.1667e+003 27.22 This problem is solved in an identical fashion to that employed in Example 27.12. For part (a), the solution is as displayed in the following plot: 6 3 0 0 10 20 30 (b) The solution for this set of equations is laid out in Sec. 28.2 (Fig. 28.9). 27.23 27.24 27.25 27.26 k=1; kmw2=[2*k,-k,-k;-k,2*k,-k;-k,-k,2*k]; [v,d]=eig(kmw2) v= 0.8034 0.1456 0.5774 -0.2757 -0.7686 0.5774 -0.5278 0.6230 0.5774 d= 3.0000 0 0 0 3.0000 0 0 0 0.0000 Therefore, the eigenvalues are 0, 3, and 3. Setting these eigenvalues equal to mω 2 , the three frequencies can be obtained. 2 mω 1 = 0 ⇒ ω 1 = 0 (Hz) 1st mode of oscillation 2 mω 2 = 0 ⇒ ω 2 = 3 (Hz) 2nd mode 2 mω 3 = 0 ⇒ ω 3 = 3 (Hz) 3rd mode 27.7 (a) The exact solution is y = Ae 5t + t 2 + 0.4t + 0.08 If the initial condition at t = 0 is 0.8, A = 0, y = t 2 + 0.4t + 0.08 Note that even though the choice of the initial condition removes the positive exponential terms, it still lurks in the background. Very tiny round off errors in the numerical solutions bring it to the fore. Hence all of the following solutions eventually diverge from the analytical solution. (b) 4th order RK. The plot shows the numerical solution (bold line) along with the exact solution (fine line). 15 10 5 0 -5 -10 0 1 2 3 4 (c) function yp=dy(t,y) yp=5*(y-t^2); >> tspan=[0,5]; >> y0=0.08; >> [t,y]=ode45('dy1',tspan,y0); (d) >> [t,y]=ode23S('dy1',tspan,y0); (e) >> [t,y]=ode23TB('dy1',tspan,y0); 30 20 10 0 -10 0 -20 1 2 3 4 -30 RK4 Analytical ODE23S ODE23TP ODE45 5