Chapter 18 VBA Program to Implement Newton interpolation

18.14 Here is a VBA program to implement Newton interpolation. It is set up to solve Example 18.5: Option Explicit Sub Newt() Dim n As Integer, i As Integer Dim yint(10) As Single, x(10) As Single, y(10) As Single Dim ea(10) As Single, xi As Single Range("a5").Select n = ActiveCell.Row Selection.End(xlDown).Select n = ActiveCell.Row - n Range("a5").Select For i = 0 To n x(i) = ActiveCell.Value ActiveCell.Offset(0, 1).Select y(i) = ActiveCell.Value ActiveCell.Offset(1, -1).Select Next i Range("e3").Select xi = ActiveCell.Value Call Newtint(x(), y(), n, xi, yint, ea) Sheets("Sheet1").Select Range("d5:f25").ClearContents Range("d5").Select For i = 0 To n ActiveCell.Value = i ActiveCell.Offset(0, 1).Select ActiveCell.Value = yint(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = ea(i) ActiveCell.Offset(1, -2).Select Next i Range("a5").Select End Sub Sub Newtint(x, y, n, xi, yint, ea) Dim i As Integer, j As Integer, order As Integer Dim fdd(10, 10) As Single, xterm As Single Dim yint2 As Single For i = 0 To n fdd(i, 0) = y(i) Next i For j = 1 To n For i = 0 To n - j fdd(i, j) = (fdd(i + 1, j - 1) - fdd(i, j - 1)) / (x(i + j) - x(i)) Next i Next j xterm = 1# yint(0) = fdd(0, 0) For order = 1 To n xterm = xterm * (xi - x(order - 1)) yint2 = yint(order - 1) + fdd(0, order) * xterm ea(order - 1) = yint2 - yint(order - 1) yint(order) = yint2 Next order End Sub 18.15 Here is the solution when the program from Prob. 18.14 is run. 18.16 See solutions for Probs. 18.1 through 18.3. 18.17 A plot of the error can easily be added to the Excel application. The following shows the solution for Prob. 18.4: The following shows the solution for Prob. 18.5: 18.18 Option Explicit Sub LagrInt() Dim n As Integer, i As Integer, order As Integer Dim x(10) As Single, y(10) As Single, xi As Single Range("a5").Select n = ActiveCell.Row Selection.End(xlDown).Select n = ActiveCell.Row - n Range("a5").Select For i = 0 To n x(i) = ActiveCell.Value ActiveCell.Offset(0, 1).Select y(i) = ActiveCell.Value ActiveCell.Offset(1, -1).Select Next i Range("e3").Select order = ActiveCell.Value ActiveCell.Offset(1, 0).Select xi = ActiveCell.Value ActiveCell.Offset(2, 0).Select ActiveCell.Value = Lagrange(x(), y(), order, xi) End Sub Function Lagrange(x, y, order, xi) Dim i As Integer, j As Integer Dim sum As Single, prod As Single sum = 0# For i = 0 To order prod = y(i) For j = 0 To order If i <> j Then prod = prod * (xi - x(j)) / (x(i) - x(j)) End If Next j sum = sum + prod Next i Lagrange = sum End Function Application to Example 18.7: 18.19 The following VBA program uses cubic interpolation for all intervals: Option Explicit Sub Newt() Dim n As Integer, i As Integer Dim yint(10) As Single, x(10) As Single, y(10) As Single Dim ea(10) As Single, xi As Single Range("a5").Select n = ActiveCell.Row Selection.End(xlDown).Select n = ActiveCell.Row - n Range("a5").Select For i = 0 To n x(i) = ActiveCell.Value ActiveCell.Offset(0, 1).Select y(i) = ActiveCell.Value ActiveCell.Offset(1, -1).Select Next i Range("e4").Select xi = ActiveCell.Value ActiveCell.Offset(2, 0).Select ActiveCell.Value = Interp(x(), y(), n, xi) Range("a5").Select End Sub Function Interp(x, y, n, xx) Dim ii As Integer If xx < x(0) Or xx > x(n) Then Interp = "out of range" Else If xx <= x(ii + 1) Then Interp = Lagrange(x, y, 0, 3, xx) ElseIf xx <= x(n - 1) Then For ii = 0 To n - 2 If xx >= x(ii) And xx <= x(ii + 1) Then Interp = Lagrange(x, y, ii - 1, 3, xx) Exit For End If Next ii Else Interp = Lagrange(x, y, n - 3, 3, xx) End If End If End Function Function Lagrange(x, y, i0, order, xi) Dim i As Integer, j As Integer Dim sum As Single, prod As Single sum = 0# For i = i0 To i0 + order prod = y(i) For j = i0 To i0 + order If i <> j Then prod = prod * (xi - x(j)) / (x(i) - x(j)) End If Next j sum = sum + prod Next i Lagrange = sum End Function Application to evaluate ln(2.5): 18.20 Sub Splines() Dim i As Integer, n As Integer Dim x(7) As Single, y(7) As Single, xu As Single, yu As Single Dim dy As Single, d2y As Single Range("a5").Select n = ActiveCell.Row Selection.End(xlDown).Select n = ActiveCell.Row - n Range("a5").Select For i = 0 To n x(i) = ActiveCell.Value ActiveCell.Offset(0, 1).Select y(i) = ActiveCell.Value ActiveCell.Offset(1, -1).Select Next i Range("e4").Select xu = ActiveCell.Value Call Spline(x(), y(), n, xu, yu, dy, d2y) ActiveCell.Offset(2, 0).Select ActiveCell.Value = yu End Sub Sub Spline(x, y, n, xu, yu, dy, d2y) Dim e(10) As Single, f(10) As Single, g(10) As Single, r(10) As Single, d2x(10) As Single Call Call Call Call Tridiag(x, y, n, e, f, g, r) Decomp(e(), f(), g(), n - 1) Substit(e(), f(), g(), r(), n - 1, d2x()) Interpol(x, y, n, d2x(), xu, yu, dy, d2y) End Sub Sub Tridiag(x, y, n, e, f, g, r) Dim i As Integer f(1) = 2 * (x(2) - x(0)) g(1) = x(2) - x(1) r(1) = 6 / (x(2) - x(1)) * (y(2) - y(1)) r(1) = r(1) + 6 / (x(1) - x(0)) * (y(0) - y(1)) For i = 2 To n - 2 e(i) = x(i) - x(i - 1) f(i) = 2 * (x(i + 1) - x(i - 1)) g(i) = x(i + 1) - x(i) r(i) = 6 / (x(i + 1) - x(i)) * (y(i + 1) - y(i)) r(i) = r(i) + 6 / (x(i) - x(i - 1)) * (y(i - 1) - y(i)) Next i e(n - 1) = x(n - 1) - x(n - 2) f(n - 1) = 2 * (x(n) - x(n - 2)) r(n - 1) = 6 / (x(n) - x(n - 1)) * (y(n) - y(n - 1)) r(n - 1) = r(n - 1) + 6 / (x(n - 1) - x(n - 2)) * (y(n - 2) - y(n - 1)) End Sub Sub Dim Dim Dim Interpol(x, y, n, d2x, xu, yu, dy, d2y) i As Integer, flag As Integer c1 As Single, c2 As Single, c3 As Single, c4 As Single t1 As Single, t2 As Single, t3 As Single, t4 As Single flag = 0 i = 1 Do If xu >= x(i - 1) And xu <= x(i) Then c1 = d2x(i - 1) / 6 / (x(i) - x(i - 1)) c2 = d2x(i) / 6 / (x(i) - x(i - 1)) c3 = y(i - 1) / (x(i) - x(i - 1)) - d2x(i - 1) * (x(i) - x(i - 1)) / 6 c4 = y(i) / (x(i) - x(i - 1)) - d2x(i) * (x(i) - x(i - 1)) / 6 t1 = c1 * (x(i) - xu) ^ 3 t2 = c2 * (xu - x(i - 1)) ^ 3 t3 = c3 * (x(i) - xu) t4 = c4 * (xu - x(i - 1)) yu = t1 + t2 + t3 + t4 t1 = -3 * c1 * (x(i) - xu) ^ 2 t2 = 3 * c2 * (xu - x(i - 1)) ^ 2 t3 = -c3 t4 = c4 dy = t1 + t2 + t3 + t4 t1 = 6 * c1 * (x(i) - xu) t2 = 6 * c2 * (xu - x(i - 1)) d2y = t1 + t2 flag = 1 Else i = i + 1 End If If i = n + 1 Or flag = 1 Then Exit Do Loop If flag = 0 Then MsgBox "outside range" End End If End Sub Sub Decomp(e, f, g, n) 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 End Sub Sub Substit(e, f, g, r, n, x) Dim k As Integer 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 18.21 The following shows the solution for Prob. 18.4: The following shows the solution for Prob. 18.5: 18.22 f1 ( x ) = f ( xo ) + f ( x1 ) − f ( xo ) ( x − xo ) x1 − xo  6.7664 − 6.5453  f1 ( x ) = 6.5453 +   ( x − 0.11144)  0.12547 − 0.11144  f1 ( x ) = 4.789107 + 15.579 x x = 0.118, f1( x) = 6.6487 s = 6.6487 kJ kg  K