Quadrature for a Mononomial

In my last post, I mentioned that my next step was to write a lift and moment routine, and that I’d been having difficulties with the interface between Windows 7 and Racket’s Scribble. Sadly, all my time (okay, all the time I’ve devoted to programming) has been devoted to Windows vs Scribble. And the loser is...(drum roll, please) ME! One of the neat things about Scribble documents are their ability to invoke live Racket commands and show their results via the interaction command. Like a good citizen in this unfriendly online world, interaction is properly suspicious of what it perceives to be foreign code. Unfortunately, my use of submodules to implement unit testing set off Racket’s paranoia. Thanks to the friendly Racket community , I eventually found call-with-trusted-sandbox-configuration, and now I can go back to scribbling and programming, at least until I fall down the next rabbit hole.

Okay, enough "mea culpa" about why it’s taken me so long to write a simple integration routine. ZFoil will use a simple enhancement of the trapezoidal rule for all its integration needs. The trapezoidal rule basically assumes that a function f is linear between points a & b, and so can estimate the integral of f by using the formula for the area of a trapezoid:

Note that this still works perfectly well when f crosses the y axis.

Now ZFoil predicts velocity pretty well, as I’ve shown. What’s more, it can predict stagnation points (where velocity touches zero, basically near the leading edge, and at the trailing edge), pretty well. A lot of the numbers we care about depend not on velocity, but of some power of velocity. Pressure (and therefore lift and moment) depends on the square of velocity, and in laminar flow, the fifth power of velocity takes on importance. Consider the plot of \(x\), \(x^2\), and a naive linear fit of \(x^2\) below:

As you can see, there’s a big difference between the area below the lines \(x^2\) and the linear fit (the green horizontal line at the top of the plot). Since ZFoil naturally spits out the sign of the velocity, we should take advantage of that information!

So, we have a function \(v(x)\). We assume \(v\) is linear between \(x_0\) and \(x_1\). We are trying to find \[\int_{x_0}^{x_1} v^n dx\] Let’s introduce the linear transformation \(T(x)=t\), such that when \(x = x_0, t =0\), and when \(x = x_1, t =1\). \[\int_{x_0}^{x_1} v^n dx = \int_0^1 v^n \frac{dx}{dt} dt = (x1 - x0)\int_0^1 v^n dt \] That last integral is easily evaluated: \[\int_0^1 v^n dt = \frac{1}{n+1}(v_1^{n+1}-v_0^{n+1})\] ...which is easily translated into a Racket function:

(define (integrate v0 v1 n)

(/(-(expt v1 (add1 n))

(expt v0 (add1 n)))

(add1 n)))

or, if you C refugees prefer, to be just a bit faster,

(define(integrate v0 v1 n)

(let[(n+1(add1 n))]

(/(-(expt v1 n+1)

(expt v0 n+1))

n+1)))

To me, that’s another plus for Racket, its flexible, and therefore expressive, variable names.

Ok, I’ve laid out my next baby step. Now I just have to take it, and incorporate unit tests, and see if I’ve truly killed the blogging bug I mentioned at the beginning of this post.