Simpson's Rules (and others)

Discussion in 'Boat Design' started by Mike D, Oct 1, 2002.

1. Joined: Sep 2002
Posts: 104
Likes: 18, Points: 0, Legacy Rep: 465

Mike DSenior Member

Simpson’s Rules are among the most useful tools that designers have in their “math toolbox”. They are the simplest to apply, the least accurate of similar methods and the easiest to screw up - when I say similar, I exclude the trapezoidal rule and counting squares.

Simpson’s First Rule (1.4.1) and Second Rule (1.3.3.1) are the second and third in a very large set of rules that have from 2 ordinates (the Trapezoidal Rule) up to as many as you want but they get more complex the more you use. They are all in the “family” of the Newton-Cotes set and Simpson’s have been in use for about 250 years.

The idea behind them all is to find the area under a closed curve; it could be a pond or a lake on a map, a boat’s waterplane or a curved sail. On a boat or a ship port and starboard are usually mirror images so it is convenient to work using only one side then double the answer.

Before we get down to it a few words of advice, when you are given or when you prepare a drawing be sure you have the right scale and the right scale-rule to make measurements. If it is a CAD drawing then do not believe the scale indicated, the output for a computer generated drawing can be on any scale, so check it before you start. It sounds obvious but I know of one case where a potential contract for 300 fast patrol boats collapsed because of a comparable basic error in the original design, there are too many similar tales to recount.

Back to Simpson.
In the past one found the area under the curve, then it became the mensuration of the curve and nowadays it has become quadrature. Therefore, if you ever go searching the net these terms may help.

Every formula in the Newton-Cotes set was derived by assuming the area could be found by one application of the rule. Look at Figs 1 and 2 below, the entire curve length is enclosed by one set of multipliers. The curves shown may seem relatively simple but they could easily have been the result of a complicated mathematical formula. Without getting into a mathematical quagmire just accept that Simpson’s First Rule is exact for a simple formula the Second Rule is exact for one a little more complex and so on. The result is, for example, that the seventh in the series handles a curve four orders higher but just see how complex it is compared to Simpson’s First;
Simpson; L/6*(1*y1 + 4*y2 + 1*y3)
Seventh: L/980*(41*y1+216*y2+27*y3+272*y4+27*y5+216*y6+41*y7)
So it explains why people would rather use repeated applications of Simpson repeated; L/18*(1*y1+4*y2+2*y3+4*y4+2*y5+4*y6+1*y7)

In these formulas the terms y1, y2 etc are the ordinates A1, A2 and A3 or B1, B2, B3 and B4 etc for the curves in the illustrations below.

If S1 (I’m tired of typing Simpson’s First Rule) is used for Fig 3 below you would get an area corresponding to a simple parabola through the three ordinates used. Instead of finding the correct area you would get the area corresponding to the dotted line. But if the curve was really a seventh order curve you would get the right answer by using the Seventh of Newton-Cotes. Using the repeating Simpson you’ll get very close. If you used more ordinates for Simpson you’d get even closer unless you requested a huge number of them in a CAD system and you might get errors due to the precision of the numbers in the calculation.

So this gives us the first clue about how many ordinates to use but first a side-bar on manual short cuts. Look at Fig 4 below, a simple curve, a “lazy” curve, not very much of it but enough to have to figure it out.
• So draw a line between the crowns, marked as AB in the figure.
• Then a perpendicular at mid-way along the base-line, marked as CD in the figure.
• Find the intersection of the line AB with the line CD.
• Measure the distance of this point from the curve.
• Above the intersection point measure 1/3 of this distance, shown with a short line in the Fig 4.
• Measure the short line position to the base, this times the base-line length is the area.
• If the curve looks like the one in Fig 4 then use the S1 method as explained later.
If you draw a line in Fig 4 through the short line and parallel to the line through the horns you will develop a parallelogram ABEF with an area and a centre of area the same as those of the curved figure ADBEF.

A standard sheer line and a standard camber line are parabolic, see Fig 5. The area is 2/3*L*H, the centre of area is 3/8*L from the left hand end and 2/5*H above the base-line. Good designers remember things like this and save time while colleagues reach for computers because they can’t see the wood for the trees.

Back to Simpson.
If we look at Fig 3 and Fig 5 (but Fig 5 is now assumed to be just another curve, we don’t know it is parabolic) we can see that there is a part where there is more curvature and this requires additional ordinates. Most people simply add another set of divisions where one part is two divisions in the 1.4.1 set. Fig 1 has only one part so they would add perhaps one or two more in Figs 3 and 5 and the multipliers would go from 1.4.1 to 1.4.2.4.1 or even 1.4.2.4.2.4.1 but this is often overkill.

Suppose Fig 5 had a small rectangle attached on the left turning it into Fig 6 then any additional ordinates along the rectangle would not help. S1 gives the correct answers on straight-line figures such as triangles and rectangles; it is only the highly curved parts that need the extra treatment. Often it is better to simply add intermediate ordinates along the stretch where they are needed which is what is shown in Fig 7.

When you add half parts, i.e. divisions, you must always keep the groups so that you have the three multipliers but because the divisions are half the basic length then the multipliers are half also. We could take, in Fig 7, ¼ ordinates at the positions shown and it would lead to more accuracy. However, Fig 5 is parabolic and it is of the same degree as the standard S1 curve so we would get the correct answer with only 3 ordinates, more ordinates only add work without benefit. But this is 20/20 hindsight, if you did not know the curve was parabolic them you should add some extra ordinates.

Look at Fig 3 again. That hump is a special case and there should be an ordinate at the peak. Two options;
• Divide the entire curve just to the right of the peak and treat it as two separate ones. Break it where the curve is still part of the down slope of the hump and not yet at the point where it is changing and forming the “hollow” on the right. The mid ordinate of this left hand part must be on or very close to hump. If it isn’t and it can’t be done you must add more ordinates. The right hand side is not too important in that it is a “gentle” and does not need as many ordinates.
• If you can’t get it to work you must simply handle the entire curve as one and add more ordinates until you have enough.
• How many is enough? Hard to tell, it depends on how abrupt that hump is and where it is. You must experiment and get a feeling for these things or you blanket it with ordinates and hope for the best. But the ultimate hump is a chine and the blanketing will miss it. The chine must be the end ordinate of one part i.e. a one of the 1.4.1, a one or a two in the 1.4.2.4.2………1 series, or a 1, 1½. ½ in a series like this 1.4.1½.2.½. This means almost certainly that you must handle the curve as two separate curves.
In almost all the net sites that describe S1 they state that an even number of divisions must be made and sometimes they talk about an even number of equally spaced divisons. This can be misleading as the previous paragraph gives a 1.4.1½.2.½ series which is 1.4.1 and ½.2.½ added on one curve. An even number of divisions but not equally spaced, the forward part is half the length of the after part.

You could break the length into say a 60/40 split if you wished. You must select one part as the main and change the multipliers in the other part in proportion. Taking the 60 as the main part then the CI (common interval) is 30% of the baseline and the multipliers are 1.4.1 in the main part and 40/60.160/60.40/60 in the other. It would be easier to use the smaller part as the main in which case the multipliers are 1½:6:1½ and 1.4.1 using a CI of 20% of the base-line length.

To be sure that you have the right multipliers and CI, add the multipliers and multiply the result by the CI as a fraction of the base length. So a simple 1.4.1 is 6*0.5 = 3. The second example in the preceding paragraph is (9+6)*.2 = 3, if you get something other than 3 it is wrong. The sum of the ordinates times their multipliers divided by the sum of the multipliers gives the average length of the ordinates. Two examples:
• Triangle: (1*0+4*0.5+1*1)/6 = 0.5 and this times the base length is the area.
• Parabola; (1*0+4*0.75+1*1)/6 = 2/3
Handling only one waterline is straight-forward but suppose you were calculating a series of waterlines for the hydrostatic curves or curves of form as they are often called.

Fig 8 shows a side profile of a curved bow with the forward three stations indicated. Even if we add half stations as shown by the intermediate lines there would still be errors in all the waterlines except the load waterline. None of them finish on the end ordinate of half parts i.e. none finish on Stn 10. We could add ¼ parts between 9 and 10; the error would reduce but it would still be present. The best way is to introduce a half station at 8 ½ so 9 is now the end of a part and the waterlines forward of 9 can all be treated as separate curves.

It isn’t possible for anyone sitting in the glory-hole to say precisely where these divisions should be made because it depends on the areas, the shapes and the degree of accuracy required. With a CAD system intermediates can be added without breaking the waterline, you must simply arrange these new intermediates so that the right hand end of a part intersects the waterline and the bow line. You will have quite a number of new ones but you’d be getting accurate answers.

On modern computers and with decent programs you can get reliable answers simply by blitzing all these difficult humps and cranks and so on. The errors become so small that they are virtually impossible to see but these tricks can’t be used in manual calculations. My opinion is that it is acceptable but poor design practice. An unusual problem arose one time and the CAD system had no direct solution module so the designer attacked it by throwing numbers at it, transferring values to another package, plotting results, and found the answer. It took three days. His boss looked at it and said it can’t be. He thought for a few minutes and sat at the designer’s computer, 20 minutes later he had the correct answer using GoalSeek in Excel.

Not all calculations integrating curves were done using Simpson’s Rules. Tchebycheff was normally used to prepare the Cross Curves of Stabilty with the integrator because it saved time, you only added the values because there are no multipliers. But the most accurate of all these rules is Gauss, the disadvantage is that the ordinates are not uniformly spaced and the odd spacings and multipliers make it very awkward to operate unless it is programmed in a CAD system.

If the system has an accurate method of determining intermediate offsets and not just a linear interpolation routine there is some justification for its use. But if the hull form has a knuckle or a chine then the curve must be segmented in the same way as Simpson. Doing it manually is not recommended due to the added risk of making an error. Following is a simple comparison of various rules to calculate the area of the tan curve up to 60 degrees using 5 ordinates, the answer correct to 5 decimals is 0.69315;
• Simpson’s First Rule (using 3 ordinates); 0.70537
• Simpson’s First Rule by its double application; 0.69452
• Newton-Cotes; 0.69379
• Tchebycheff (or Chebyshev); 0.69289
• Gauss; 0.69317

Simpson’s Second Rule is not too much used, it can be handy for quick manual jobs but S1 suits the standard ordinate spacing much better. S2 is slightly more accurate but with special naval architectural programs available as well as the excellent spreadsheets it is effortless to simply increase the number of ordinates. Given a choice it is wiser to use S1 and thoroughly learn its pitfalls which exist in varying degrees in all these methods of integration.

Suppose that a double application of S1 is considered “normal” with its error of approximately 0.2%. How does it compare with short waterlines as shown in Fig 7? It’s difficult to be sure because a lot depends on how many ordinates were used. The guidance values that I’ll demonstrate seem far too much but look at it this way.

Imagine a drawing of waterlines in the forebody with stations from 5 to 9 in halves then from 9 to 10 in quarters. If the waterline stops on Stn 9 ½ you would get the correct answer when the line is curved. But when it stops elsewhere there is an error because S1 assumes that the curve extends over the three equi-spaced ordinates so it assumes a curve through the first declared points as given and then to zero at the others. Therefore, the curve it uses does not resemble the real one so the errors as caused. There are four charts shown below that illustrate this, each for a different waterplane area coefficient but based on the distance from Station 5 to where the curve meets the centre-line.

The x-axis of the charts is where the curve meets the centre-line and is shown only from Station 9 to Station 10. The lines are denoted as ½’s, 1/4’s etc meaning the divisions so that the multipliers are 1.4.1 for ½’s, 1.4.2.4.1 for ¼’s, 1.4.2.4.2.4.2.4.1 for 1/8’s and 1.4.11/2.2.1.2.1/2 for ¼’s and 1/8’s. The last one is in quarters from Stn 9 to 9.5 and in eighths from Stn 9.5 to 10. The first three charts are almost look-alikes and the growth in the values is only a result of the increasing Cw, but you can see the similarity of the results until we look at the last chart.

With a Cw of 1 it is a rectangular shape and it could be the afterbody with a wide transom. Another factor kicks in because the half-breadth at the transom is not zero. Say the after body used S1 with the stations in halves, so that would be 0, ½, 1, 1½ etc. Then put the transom on Stn 1, but the multiplier on Stn 1 is doubled because S1 assumes one curve fro 0 to 1 and another from 1 to 2 and so on. Remember the check I showed earlier about the sum of the multipliers? If we add them from 1 to 5 we get 24 but as the multipler at Stn 1 is doubled we actually used 25 so the error is 1/24 or 4.167% and that’s what is plotted.

So you see how easy it is to make a mistake and I’m quite sure you did not realize quite how big these errors can be. Once again, if you need help just ask!

Michael

Attached Files:

• simpson.gif
File size:
164.5 KB
Views:
5,091
2 people like this.
2. Joined: Sep 2001
Posts: 59
Likes: 8, Points: 8, Legacy Rep: 82
Location: Rhode Island

Steve HollisterJunior Member

Another way to look at integration is as a combination of curve fitting followed by integration. For calculations like a hull section's area (not sectional area curve) you often have non-simple station shapes that have either two values for each half breadth or two values for each height value. These multi-valued function shapes cannot be integrated using Simpson-type rules. The old Navy SCHP hull calculation program forced one into geometry contortions because of this problem.

The solution, which has been around for ages, is to integrate using line integrals which does just what the planimeter does. It's kind of amazing that just by tracing around a closed loop you can determine the area of the closed loop. In fact, using Green's Theorem, you can calculate areas, centroids, and more just by "marching" along a closed loop polyline. The program that does this is pretty simple, but it is still not as good as Simpson's rule if you have to do a quick and dirty integration in a spreadsheet. Then again, how difficult is it to type in the coordinates of the curve that your are integrating into a CAD program? I believe that there are add-ins to AutoCAD (and perhaps other programs) to do these calculations.

The nice part of the line integral approach in a CAD program is that the curve-fitting is separated from the integration. In other words, if the closed curve looks "good" in the CAD program, then you know that the results are accurate. That is why we used this line integral approach in our calculation software, including the longitudinal integration of the sectional area curve. If the displayed curve looks good, then the integration is good.

Steve Hollister

1 person likes this.
3. Joined: Sep 2002
Posts: 104
Likes: 18, Points: 0, Legacy Rep: 465

Mike DSenior Member

Steve

I grinned at visions of contortions due to simplicity. It's about 40 years since I used Green and it was only to try it out - long-hand. I wasn't impressed because I was able to work out the forms you describe faster just using a slide-rule and a scratch pad. Now that's giving away my age, who on this forum has ever used one I wonder.

All years before computers, of course, but nowadays with the proliferation of PC's and software costing a fraction of the early issues it isn't surprising that Green is being adopted.

You're right about planimeters, I was amazed when I first used one and when I was allowed to use the Amsler Integrator!! First and second moments, by God.

I liked your reference to down and dirty spreadsheets, I like the ability to see what's going on and to be able to do simple GoalSeeking and Solver problems.

The principal reason for my posting was to pass on some of my experience with Simpson's Rules and to point out some of the strange things that happen, invariably unexpectedly and at the most awkward moments.

It was refreshing to read of the attitude you have adopted in your software development - looking good being good. There are too many naive souls still believing that if it comes from a computer it must be good. I wrote a short program some years ago with a serious, deliberate mistake in it. No one found it because the ouput sheet was beautifully formatted, a real delight. But now I'm on my old hobby-horse - Form vs Function.

Michael

4. Joined: Apr 2002
Posts: 204
Likes: 1, Points: 0, Legacy Rep: 16
Location: Australia

BrettMSenior Member

Steve,

The Autocad Function is under region/mass properties (command=massprop) It is standard even with AutoCad LT.

Brett

5. Joined: Sep 2001
Posts: 59
Likes: 8, Points: 8, Legacy Rep: 82
Location: Rhode Island

Steve HollisterJunior Member

Mike D,

For me, it all started when I was at U.of Michigan in the mid 70's (I don't feel that old!) tearing apart the Navy's SHCP program (Ship's Hull Characteristics Program - FORTRAN spaghetti code with huge EQUIVALENCE blocks. How many remember those things?) which used an odd, sequential, 3 point parabolic fit of station offsets to do the station integration vertically. Someone had written a Calcomp plotter routine to plot the stations to see if the input data looked OK. However, the plot program used a spline function to fit the offsets. Therefore, the plotted stations could look great, but the 3 point parabolic fit and the integration could be quite wrong! Oops!

By the way, we underclassmen weren't allowed to touch the integrator. (I don't remember which model it was.) I do remember wondering how on earth it could possibly calculate areas and moments just by following a closed path. I started college using slide rulers and finished using computers. (I even had a circular slide ruler once!)

Steve

6. Joined: Sep 2002
Posts: 104
Likes: 18, Points: 0, Legacy Rep: 465

Mike DSenior Member

Steve

I think the college authorities were quite right in keeping you young rowdies away from valuable artifices.

But I can now admit, in deep humility, that I was once allowed to touch our sister company's integraph but not use it. Three machines have made me gape in awe - Babbage's Analytical Engine, Kelvin's Tide Predictor and the Integraph and I got to touch one, I can still remember the shiver as I watched someone use it. WOW!!!!

Michael

7. Joined: Oct 2001
Posts: 3,590
Likes: 130, Points: 0, Legacy Rep: 2369
Location: Australia

WillallisonSenior Member

Sometimes, you guys amaze me - wish I knew what the bloody hell you are talking about!!

8. Joined: Dec 2002
Posts: 8
Likes: 0, Points: 0, Legacy Rep: 10
Location: Oahu, Hawaii

BrianJunior Member

SHCP & Computer Output

I appreciate the insightful comment made by Steve regarding code inconsistencies in SHCP (i. e. the section plotting function is not a true representation of the calculated function).

The comments by Mike that computer output is often accepted as being correct even if it is wrong. This demonstrates a weaknesses in human nature.

I think it is really important that people understand how their computer software is calculating things. They should know how to calculate the results by hand (with the aid of a calculator, planimeter, integrator, spreadsheet or CAD program) or use a different program to verify the results (at least when starting to use a new program). They should be at least be aware of simple rules of thumb methods to quickly check computer output.

I think programs are great because they can do tons of calculations without getting tired. But if one cannot understand the process and check the program's output it is dangerous to use the program.

9. Joined: Aug 2002
Posts: 107
Likes: 2, Points: 0, Legacy Rep: 64
Location: Crystal River, FL USA

TimmSenior Member

I second that. Most of the programs I use are spreadsheets that are written from the same formulas I used before I got a computer. I always do a complete calculation by hand (so to speak) and then check it with any new spreadsheet to make sure I got all the formulas in all the little cells correctly.

10. jsalGuest

Hello my name is John and I am a student of Naval Academy of Greece. As I can read above you are very familiar with Simpson's Rules. I have a report to prepare about Simpson's Rules and its use on the calculation of the area and the centre of weight of ships. Though I succeed in finding information about the calculation of the waterplane area I fear that I am unable to the same for the centre of weight. I wonder if you are so kind to share some knowledge with me providing me with some information. You can mail me at the salalalas@hotmail.com Thank you for your time. With all my respect John.

11. Joined: Feb 2004
Posts: 223
Likes: 1, Points: 0, Legacy Rep: 15
Location: Connecticut

As I recall, the error term for S1 involves the square of then interval. Thus halving the interval will reduce the error by 75%. (This would be error for smooth but non-parabolic curves.)

When programming this sort of thing, it is simpler to divide the shape into triangles, but care should be taken since the usual formulae are subject to calculation error for long, skinny triangles. Many formulas familar from plane geometry are inferior to simple expressions derived from analytic (Cartesian) geometry (e.g. for the centroid of a triangle).

How does one decide when the internal accuracy of the drawing is better than the difference between the drawing and the actual boat?

12. Joined: Feb 2002
Posts: 511
Likes: 27, Points: 28, Legacy Rep: 394
Location: Denmark

sorenfdkYacht Designer

The attachments to Mikes original post seem to be missing. I'd really like to see them, so if anybody has them, please post them here (again!)

Søren

13. Joined: Jun 2001
Posts: 1,368
Likes: 71, Points: 58, Legacy Rep: 923
Location: Great Lakes

JeffModerator

The figures are still attached to the original post above (click on the thumbnail). I don't remember there being more attachments, but if you remember seeing more let me know and I'll have a look.

14. Joined: Feb 2002
Posts: 511
Likes: 27, Points: 28, Legacy Rep: 394
Location: Denmark

sorenfdkYacht Designer

I can't see any thumpnail to click, but anyway I got them now.

Thanks!

Søren

Forum posts represent the experience, opinion, and view of individual users. Boat Design Net does not necessarily endorse nor share the view of each individual post.
When making potentially dangerous or financial decisions, always employ and consult appropriate professionals. Your circumstances or experience may be different.