Hi. This should do it. This is something I haven't yet finished up. You
might also do a Find for "HOWDY" to go beyond trapezoidal.
HTH
Dave Braden
A Worksheet-based Solution
Suppose your x-values are in a column-range you named Xpts, and the
corresponding y-values are in a columns-range Ypts. What is the area
under the curve? To keep the following worksheet formula short, choose
Insert/Name/Define, and set "Pnls" to be =COUNT(Xpts)-1. This is the
number of panels formed by successive pairs of (x1,x2) points.
Alternatively, you can enter the formula into a cell, and then
Insert/Name/Define that cell to be Pnls.
Next, enter this function into a convenient cell:
=SUMPRODUCT((OFFSET(Xpts,1,0,Pnls)-OFFSET(Xpts,0,0,Pnls)),(OFFSET(Ypts,1,0,Pnls)+OFFSET(Ypts,0,0,Pnls)))*0.5
The values in Xpts must be nondecreasing (e.g., 1.2, 2.3, 2.3, 3.0, 4.0)
or all non-increasing (e.g., 5.2, 4.0, 4.0, 3.1, 2.5). You can build in
some error-checking for this by modifying the above formula to
=SUMPRODUCT((OFFSET(Xpts,1,0,Pnls)-OFFSET(Xpts,0,0,Pnls)),(OFFSET(Ypts,1,0,Pnls)+OFFSET(Ypts,0,0,Pnls)))*0.5*IsWeakMonotone
where "IsWeakMonotone" is the name of a cell containing the
array-entered formula
=IF(OR(AND(OFFSET(Xpts,1,0,Pnls)>=OFFSET(Xpts,0,0,Pnls)),AND(OFFSET(Xpts,1,0,Pnls)<=OFFSET(Xpts,0,0,Pnls))),1,NA())
Remember, to enter an array formula, hold ctrl-shift when pressing Enter.
This is a (theoretically) exact solution to the problem. It sums up the
areas of the individual trapezoids formed by the panels. To get a
graphic feel for this, and where we're headed below, try it!
If Xpts forms a *row* range, adjust the corresponding OFFSETs above for
Xpts to OFFSET(Xpts,0,1,0,Pnls) and OFFSET(Xpts,0,0,0,Pnls). Do
similarly for Ypts as needed. If one is a row and the other a column,
use Excel's TRANSPOSE function (see Excel's Help if you are unfamiliar
with this); I'd recommend using it on the row-oriented range in the
Named Definition.
A VBA Solution
If you want a more flexible solution, you can use the following VBA
(Visual Basic for Applications) code. It checks for weak monotonicity of
Xpts, and other potential problems. This code has been tested for Excel
versions 9 and later; earlier versions of Excel may require some
editing. To use, enter into a convenient cell
=TrapezIntegrate(Xpts,Ypts) where each of the arguments is either a
single-row or single-column range of cells that evaluate to numbers
Optional arguments let you determine the area under part of the curve.
Option Explicit
'////////////////////////////////////////////////////////////////////////////////
' TrapezIntegrate
' For given arrays of (x,y) pts, determine integral,
' assuming pts connected by straight lines, and that
' the x-values are nondecreasing.
'
' David J. Braden
' 2002 March 13
'////////////////////////////////////////////////////////////////////////////////
Function TrapezIntegrate(Xin As Variant, Yin As Variant, _
Optional LowX As Double = -1E+307, Optional HighX As Double =
1E+307 _
) As Variant
Dim lNumPts As Long, l As Long
Dim lLoPart As Long, lHiPart As Long
Dim bNegInt As Boolean
Dim dDelta As Double, dSum As Double
Dim Xpts, Ypts 'variants
On Error Resume Next
If Not MakeVect(Xin, Xpts) Then
TrapezIntegrate = Xpts 'return w/ error
Exit Function
End If
lNumPts = UBound(Xpts)
If Not MakeVect(Yin, Ypts) Then
TrapezIntegrate = Ypts 'return w/ error
Exit Function
ElseIf (lNumPts <> UBound(Ypts)) Or (lNumPts < 2) Then
TrapezIntegrate = CVErr(xlErrValue)
Exit Function
End If
If LowX > HighX Then 'switch them
dDelta = HighX: HighX = LowX: LowX = dDelta: bNegInt = True
End If
If HighX < Xpts(1) Or LowX > Xpts(lNumPts) Then
TrapezIntegrate = CVErr(xlErrNum)
Exit Function
End If
If LowX < Xpts(1) Then
LowX = Xpts(1)
lLoPart = 1
Else
lLoPart = Application.WorksheetFunction.Match(LowX, Xpts, 1)
End If
If HighX > Xpts(lNumPts) Then
HighX = Xpts(lNumPts)
lHiPart = lNumPts - 1
Else
lHiPart = Application.WorksheetFunction.Match(HighX, Xpts, 1)
If lHiPart = lNumPts Then lHiPart = lHiPart - 1
End If
'check that Xpta and Ypts numeric outside of range to be integrated
'the interior pots will be checked during integration
With Application.WorksheetFunction
For l = 1 To lLoPart
If Not (.IsNumber(Xpts(l)) And .IsNumber(Ypts(l))) Then
TrapezIntegrate = CVErr(xlErrValue)
Exit Function
End If
Next
For l = lHiPart + 1 To lNumPts
If Not (.IsNumber(Xpts(l)) And .IsNumber(Ypts(l))) Then
TrapezIntegrate = CVErr(xlErrValue)
Exit Function
End If
Next
End With
If lLoPart = lHiPart Then 'points are in the same partition
dSum = (LinTerp(LowX, Xpts(lLoPart), Xpts(lLoPart + 1),
Ypts(lLoPart), Ypts(lLoPart + 1)) _
+ LinTerp(HighX, Xpts(lLoPart), Xpts(lLoPart + 1),
Ypts(lLoPart), Ypts(lLoPart + 1))) _
* (HighX - LowX)
Else
dSum = (LinTerp(LowX, Xpts(lLoPart), Xpts(lLoPart + 1),
Ypts(lLoPart), Ypts(lLoPart + 1)) _
+ Ypts(lLoPart + 1)) * (Xpts(lLoPart + 1) - LowX) _
+ (LinTerp(HighX, Xpts(lHiPart), Xpts(lHiPart + 1),
Ypts(lHiPart), Ypts(lHiPart + 1)) _
+ Ypts(lHiPart)) * (HighX - Xpts(lHiPart))
For l = lLoPart + 1 To lHiPart - 1
dDelta = Xpts(l + 1) - Xpts(l)
dSum = dSum + dDelta * (Ypts(l + 1) + Ypts(l))
If Err.Number <> 0 Then 'one of the (X,Y) not numeric
TrapezIntegrate = CVErr(xlErrValue)
Exit Function
ElseIf Sgn(dDelta) < 0 Then
TrapezIntegrate = CVErr(xlErrNum)
Exit Function
End If
Next
End If
If bNegInt Then
TrapezIntegrate = CVar(-dSum * 0.5)
Else
TrapezIntegrate = CVar(dSum * 0.5)
End If
End Function
'////////////////////////////////////////////////////////////////////////////////
Private Function LinTerp(x, x1, x2, y1, y2)
On Error Resume Next
If x1 = x2 Then
If x = x1 Then LinTerp = y2 Else LinTerp = CVErr(xlErrNum)
Else
LinTerp = y1 + (y2 - y1) * (x - x1) / (x2 - x1)
End If
If Err.Number <> 0 Then LinTerp = CVErr(xlErrValue)
End Function
'////////////////////////////////////////////////////////////////////////////////
' Makevect takes an array or range, and orients it as a row vector
' Return error if vSrc has more than 1 row and more than 1 col.
'////////////////////////////////////////////////////////////////////////////////
Public Function MakeVect(vSrc As Variant, vx As Variant) As Boolean
Dim lCol As Long, lRow As Long
If IsError(vSrc) Then vx = vSrc: Exit Function
If TypeName(vSrc) = "Range" Then
If vSrc.Areas.Count > 1 Then vx = CVErr(xlErrValue): Exit Function
vx = vSrc.Value
Else
vx = vSrc
End If
On Error Resume Next
lRow = UBound(vx, 1): lCol = UBound(vx, 2)
If lCol = 1 Then 'already have a column
MakeVect = True
vx = Application.WorksheetFunction.Transpose(vx)
ElseIf lCol = 0 Then 'have a scalar or row
MakeVect = True
If lRow = 0 Then ReDim vx(1 To 1): vx(1) = vSrc
Else 'have at least 2 columns
If lRow > 1 Then
vx = CVErr(xlErrValue)
Else
With Application.WorksheetFunction
vx = .Transpose(vx)
vx = .Transpose(vx)
End With
MakeVect = True
End If
End If
End Function
B) The Area with Smoothing
Microsoft Excel apparently uses Bezier smoothing for its charts. The
downside of this, as far as the topic at hand goes, is that even if the
support of the function (the "X points") is monotone (i.e.,
nondecreasing or nonincreasing), the curve might bend back on itself,
rendering the question "what is the area under the curve" meaningless
(or at best ambiguous). To see this, enter into A1:A5 the values 1, 2,
3, 3.1, 4.5, and into B1:B5 enter 1.5, 2, 1, 3, 2.5. Select A1:B5 and
create a scatter (XY) chart. Notice how the curve "bends back" between 3
and 3.1.
However, the Bezier smooth is but one of many. There are two in
particular that warrant attention: the cubic spline, and the (double)
parabolic.
1) Cubic Spline
The cubic spline is very smooth; it is far better behaved than the
Bezier in maintaining monotonicity, yet tends, by comparison, to
"overshoot" under when, say, the data pairs take a sudden change, as in
Xpts = {1,2,2.1,3}, Ypts = {2,2,5,3}
2. I Know the Function I am Plotting
Easily the fastest way to do this within Excel is with worksheet
functions, but see below for VBA solutions. Reference material follows.
HOWDY
A) Worksheet Solutions
Here's a way to set up your worksheet to get three popular solutions,
all of which fall in the Newton-Cotes family. They are the Trapezoidal
Rule, Simpson's Rule, and Weddle's Rule.
To use Excel for evaluating the integral of, say, 2+3*(Ln(x))^0.6:
¥ Enter some labels:
In cell A1, enter "X_1"
A2 "X_n"
A3 "NbrPanels"
¥ Set some values:
In cell B1, enter 1
B2 2.5
B3 600
¥ Define some names:
Select A1:B3, then choose Insert/Name/Create. Make sure that only the
"Left Column" box is selected. If it isn't, you might have entered text
instead of numbers in the right column, or mis-selected the range. Press OK.
¥ Choose Insert/Names/Define
Enter each of the following names and their definitions, pressing Add
with each entry (you can copy and paste these):
EPanels =6*ROUNDUP(NbrPanels/6,0)
delta =(X_n-X_1)/EPanels
Steps =ROW(INDIRECT("1:"&EPanels+1))-1
EvalPts =X_1+delta*Steps
WeddleWts
=IF(MOD(Steps,EPanels)=0,1,IF(MOD(Steps,6)=0,2,IF(MOD(Steps,3)=0,6,IF(MOD(Steps,2)=0,1,5))))*(0.3*Delta)
EPanels takes NbrPanels and rounds it up to the next multiple of 6. This
isn't neccessary for Simpson's rule, which would need =EVEN(NbrPanels)
instead, or the Trapezoidal, which needs no adjustment. It *is*
required, though, for the more sophisticated Weddle's Rule.
(optional) If you are interested in a trapezoidal approximation, define
TrapWts =IF(MOD(Steps,EPanels)=0,0.5*Delta,Delta)
For Simpsons's Rule, define
SimpsonWts
=IF(MOD(Steps,EPanels)=0,1,IF(MOD(Steps,2)=1,4,2))*(Delta/3)
¥ Close the Define Names box, and in, say, cell D1, array-enter
=SUM(WeddleWts*(2+3*LN(EvalPts)^0.6))
That is, type in the function, and hold Ctl-Shift when pressing Enter.
In general, ctrl-shift-enter
=SUM(MyWts*f(EvalPts))
where f() is a legitimate Excel expression that yields a scalar numeric
value, and MyWts is any one of TrapWts, SimpsonWts or WeddleWts.
For the sample problem, Mathematica gives, at a higher precision than
Excel works with, 5.94249912704062.
Notes:
¥ NbrPanels should be at least 6 for this implementation. You typically
achieve greater accuracy as NbrPanels increases.
¥ You (usually) do better, for given NbrPanels, with SimpsonWts than
TrapWts, and better yet with WeddleWts. I don't know of a function, *as
graphed by Excel*, which isn't best handled by WeddleWts.
¥ If you have "jumps" in your function, break it up at the points
where those occur, and add the pieces.
¥ To keep keep your worksheet so that the calculations update with
changes in X_1, X_n or NbrPanels, enter =Delta in some cell (say, $B$4),
and be sure that Calculation isn't set to Manula in Excel's Preferences.
¥ Using a high value for NbrPanels (greater than 6^4 = 1296 or so) might
noticably slow down your workbook's performance under certain
circumstances. Chances are, especially with WeddleWts, you will need far
fewer points to get your desired accuracy.
B) VBA Solutions
While VBA solutions are perhaps slower than worksheet-based ones, they
offer more flexibility.
If you want a one-time solution, then the following is quick and easy.
Paste the code below into a VBA module, and run CalcInt (note-this uses
a slightly modified version of Karl's post):
'//////////////////
Option Explicit
Sub CalcInt()
Dim x As Double
x = Romberg(0, 2.5, 13)
Debug.Print " x= " & x
End Sub
' Hard-wire your function here
Function fct(x As Double) As Double
fct = 2 + 3 * (Log(x) ^ 0.6)
End Function
' The function Romberg calculates the value of the integral of
' the function fct() between the limits lgr and rgr.
' The order of the error term is 2*Order+2, Order>=0.
' Karl M. Syring, 1/6/01 (
[email protected])
Function Romberg(LowX As Double, HiX As Double, Order As Long) As Double
ReDim t(1 To Order + 1) As Double
Dim dSup As Double, dU As Double, dM As Double
Dim f As Long, h As Long, j As Long, n As Long
dSup = HiX - LowX
t(1) = (fct(LowX) + fct(HiX)) * 0.5
n = 2
For h = 1 To Order
dU = 0#
dM = dSup / n
For j = 1 To n - 1 Step 2
dU = dU + fct(LowX + j * dM)
Next j
t(h + 1) = dU / n + t(h) * 0.5
f = 1
For j = h To 1 Step -1
f = 4 * f
t(j) = t(j + 1) + (t(j + 1) - t(j)) / (f - 1)
Next j
n = 2 * n
Next h
Romberg = t(1) * dSup
End Function
'//////////////////
Running the above (where the third argument of Romberg = 13) will give,
for the particular example specified in "fct", a result comparable to
the worksheet-based approaches using NbrPanels =2*6^5; the relative
error is about 2E-08 among the methods.
For more general situations (e.g., not needing to "hard-code" the
function definition), and more advanced algorithms, contact the author
of this note.
3. References
Press, et al., Numerical Recipes in C, available on-line at
http://lib-www.lanl.gov/numerical/bookcpdf.html
Eric Weisstein's "World of Mathematics" web page
