Chapter 7 VBA Program to Implemen tthe Müller Algorithm

7.6 Errata in Fig. 7.4; 6th line from the bottom of the algorithm: the > should be changed to >= IF (dxr < eps*xr OR iter >= maxit) EXIT Here is a VBA program to implement the Müller algorithm and solve Example 7.2. Option Explicit Sub TestMull() Dim maxit As Integer, iter As Integer Dim h As Single, xr As Single, eps As Single h = 0.1 xr = 5 eps = 0.001 maxit = 20 Call Muller(xr, h, eps, maxit, iter) MsgBox "root = " & xr MsgBox "Iterations: " & iter End Sub Sub Muller(xr, h, eps, maxit, iter) Dim Dim Dim Dim x0 As Single, x1 As Single, x2 As Single h0 As Single, h1 As Single, d0 As Single, d1 As Single a As Single, b As Single, c As Single den As Single, rad As Single, dxr As Single x2 = xr x1 = xr + h * xr x0 = xr - h * xr Do iter = iter + 1 h0 = x1 - x0 h1 = x2 - x1 d0 = (f(x1) - f(x0)) / h0 d1 = (f(x2) - f(x1)) / h1 a = (d1 - d0) / (h1 + h0) b = a * h1 + d1 c = f(x2) rad = Sqr(b * b - 4 * a * c) If Abs(b + rad) > Abs(b - rad) Then den = b + rad Else den = b - rad End If dxr = -2 * c / den xr = x2 + dxr If Abs(dxr) < eps * xr Or iter >= maxit Then Exit Do x0 = x1 x1 = x2 x2 = xr Loop End Sub Function f(x) f = x ^ 3 - 13 * x - 12 End Function 7.7 The plot suggests a root at 1 6 4 2 0 -1 -2 0 1 2 -4 -6 Using an initial guess of 1.5 with h = 0.1 and eps = 0.001 yields the correct result of 1 in 4 iterations. 7.8 Here is a VBA program to implement the Bairstow algorithm and solve Example 7.3. Option Explicit Sub PolyRoot() Dim n As Integer, maxit As Integer, ier As Integer, i As Integer Dim a(10) As Single, re(10) As Single, im(10) As Single Dim r As Single, s As Single, es As Single n = 5 a(0) = 1.25: a(1) = -3.875: a(2) = 2.125: a(3) = 2.75: a(4) = -3.5: a(5) = 1 maxit = 20 es = 0.01 r = -1 s = -1 Call Bairstow(a(), n, es, r, s, maxit, re(), im(), ier) For i = 1 To n If im(i) >= 0 Then MsgBox re(i) & " + " & im(i) & "i" Else MsgBox re(i) & " - " & Abs(im(i)) & "i" End If Next i End Sub Sub Bairstow(a, nn, es, rr, ss, maxit, re, im, ier) Dim Dim Dim Dim Dim iter As Integer, n As Integer, i As Integer r As Single, s As Single, ea1 As Single, ea2 As Single det As Single, dr As Single, ds As Single r1 As Single, i1 As Single, r2 As Single, i2 As Single b(10) As Single, c(10) As Single r = rr s = ss n = nn ier = 0 ea1 = 1 ea2 = 1 Do If n < 3 Or iter >= maxit iter = 0 Do iter = iter + 1 b(n) = a(n) b(n - 1) = a(n - 1) + r c(n) = b(n) c(n - 1) = b(n - 1) + r For i = n - 2 To 0 Step Then Exit Do * b(n) * c(n) -1 b(i) = a(i) + r * b(i + 1) + s * b(i + 2) c(i) = b(i) + r * c(i + 1) + s * c(i + 2) Next i det = c(2) * c(2) - c(3) * c(1) If det <> 0 Then dr = (-b(1) * c(2) + b(0) * c(3)) / det ds = (-b(0) * c(2) + b(1) * c(1)) / det r = r + dr s = s + ds If r <> 0 Then ea1 = Abs(dr / r) * 100 If s <> 0 Then ea2 = Abs(ds / s) * 100 Else r = r + 1 s = s + 1 iter = 0 End If If ea1 <= es And ea2 <= es Or iter >= maxit Then Exit Do Loop Call Quadroot(r, s, r1, i1, r2, i2) re(n) = r1 im(n) = i1 re(n - 1) = r2 im(n - 1) = i2 n = n - 2 For i = 0 To n a(i) = b(i + 2) Next i Loop If iter < maxit Then If n = 2 Then r = -a(1) / a(2) s = -a(0) / a(2) Call Quadroot(r, s, r1, i1, r2, i2) re(n) = r1 im(n) = i1 re(n - 1) = r2 im(n - 1) = i2 Else re(n) = -a(0) / a(1) im(n) = 0 End If Else ier = 1 End If End Sub Sub Quadroot(r, s, r1, i1, r2, i2) Dim disc disc = r ^ 2 + 4 * s If disc > 0 Then r1 = (r + Sqr(disc)) / 2 r2 = (r - Sqr(disc)) / 2 i1 = 0 i2 = 0 Else r1 = r / 2 r2 = r1 i1 = Sqr(Abs(disc)) / 2 i2 = -i1 End If End Sub 7.9 See solutions to Prob. 7.5 7.10 The goal seek set up is The result is 7.11 The goal seek set up is shown below. Notice that we have named the cells containing the parameter values with the labels in column A. The result is 63.649 kg as shown here: 7.12 The Solver set up is shown below using initial guesses of x = y = 1. Notice that we have rearranged the two functions so that the correct values will drive them both to zero. We then drive the sum of their squared values to zero by varying x and y. This is done because a straight sum would be zero if f1(x,y) = -f2(x,y). The result is 7.13 A plot of the functions indicates two real roots at about (−1.5, 1.5) and (−1.5, 1.5). 4 3 2 1 0 -3 -2 -1 0 -1 -2 -3 1 2 3 The Solver set up is shown below using initial guesses of (−1.5, 1.5). Notice that we have rearranged the two functions so that the correct values will drive them both to zero. We then drive the sum of their squared values to zero by varying x and y. This is done because a straight sum would be zero if f1(x,y) = -f2(x,y). The result is For guesses of (1.5, 1.5) the result is (1.6004829, 1.561556). >> roots (d) with the expected result that the remaining roots of the original polynomial are found ans = 8.0000 -4.0000 1.0000 We can now multiply d by b to come up with the original polynomial, >> conv(d,b) ans = 1 -9 -20 204 208 -384 Finally, we can determine all the roots of the original polynomial by >> r=roots(a) r = 8.0000 6.0000 -4.0000 -2.0000 1.0000 7.15 p=[0.7 -4 6.2 -2]; roots(p) ans = 3.2786 2.0000 0.4357 p=[-3.704 16.3 -21.97 9.34]; roots(p) ans = 2.2947 1.1525 0.9535 p=[1 -2 6 -2 5]; roots(p) ans = 1.0000 1.0000 -0.0000 -0.0000 + + - 2.0000i 2.0000i 1.0000i 1.0000i 7.16 Here is a program written in Compaq Visual Fortran 90, PROGRAM Root Use IMSL !This establishes the link to the IMSL libraries Implicit None !forces declaration of all variables Integer::nroot Parameter(nroot=1) Integer::itmax=50 Real::errabs=0.,errrel=1.E-5,eps=0.,eta=0. Real::f,x0(nroot) ,x(nroot) External f Integer::info(nroot) Print *, "Enter initial guess" Read *, x0 Call ZReal(f,errabs,errrel,eps,eta,nroot,itmax,x0,x,info) Print *, "root = ", x Print *, "iterations = ", info End Function f(x) Implicit None Real::f,x f = x**3-x**2+2*x-2 End The output for Prob. 7.4a would look like Enter initial guess .5 root = 1.000000 iterations = 7 Press any key to continue 7.17 ho = 0.55 – 0.53 = 0.02 h1 = 0.54 – 0.55 = -0.01 δo = 58 – 19 = 1950 0.55 – 0.53 δ1 = 44 – 58 = 1400 0.54 – 0.55 a= δ1– δo = -55000 h1 + ho b = a h1 + δ1 = 1950 c = 44 b 2 − 4ac = 3671.85 t o = 0.54 + − 2( 44) = 0.524 s 1950 + 3671.85 Therefore, the pressure was zero at 0.524 seconds. 7.18 I) Graphically: EDU»C=[1 3.6 0 -36.4];roots(C) ans = -3.0262+ 2.3843i -3.0262- 2.3843i 2.4524 The answer is 2.4524 considering it is the only real root. II) Using the Roots Function: EDU» x=-1:0.001:2.5;f=x.^3+3.6.*x.^2-36.4;plot(x,f);grid;zoom By zooming in the plot at the desired location, we get the same answer of 2.4524. x 10 -4 4 2 0 -2 -4 -6 2.4523 2.45232.4524 2.4524 2.4524 2.4524 2.45242.4524 2.4524 2.4524 2.4524 7.19 Excel Solver Solution: The 3 functions can be set up as roots problems: f1(a , u , v) = a 2 − u 2 + 3v 2 = 0 f2 (a, u , v) = u + v − 2 = 0 f 3 ( a, u , v) = a 2 − 2 a − u = 0 Symbolic Manipulator Solution: >>syms a u v >>S=solve(u^2-3*v^2-a^2,u+v-2,a^2-2*a-u) >>double (S.a) ans = 2.9270 + 0.3050i 2.9270 – 0.3050i -0.5190 -1.3350 >>double (S.u) ans = 2.6203 + 1.1753i 2.6203 – 1.1753i 1.3073 4.4522 >>double (S.v) ans = -0.6203 + 1.1753i -0.6203 – 1.1753i 0.6297 -2.4522 Therefore, we see that the two real-valued solutions for a, u, and v are (-0.5190,1.3073,0.6927) and (-1.3350,4.4522,-2.4522). 7.20 The roots of the numerator are: s = -2, s = -3, and s = -4. The roots of the denominator are: s = -1, s = -3, s = -5, and s = -6. G (s ) = ( s + 2)(s + 3)(s + 4) ( s + 1)(s + 3)(s + 5)(s + 6)