How to represent [airfoil] coordinates as a polynomial

I am creating some code that allows the user to input a file containing airfoil coordinates ( :,1) = x coordinates, ( :,2) = y coordinates. Below are the coordinates for the upper surface of the airfoil (Leading edge to Trailing edge)

0	0
0.00369000000000000	0.0134500000000000
0.0170200000000000	0.0280700000000000
0.0336900000000000	0.0390000000000000
0.0515300000000000	0.0479600000000000
0.0699500000000000	0.0556700000000000
0.107870000000000	0.0682100000000000
0.127140000000000	0.0734700000000000
0.146570000000000	0.0781200000000000
0.166120000000000	0.0822000000000000
0.185780000000000	0.0857400000000000
0.205520000000000	0.0888000000000000
0.225320000000000	0.0914200000000000
0.245170000000000	0.0936200000000000
0.265070000000000	0.0954000000000000
0.284990000000000	0.0968000000000000
0.304940000000000	0.0978000000000000
0.324910000000000	0.0984500000000000
0.344880000000000	0.0987500000000000
0.364850000000000	0.0987300000000000
0.384830000000000	0.0984200000000000
0.404790000000000	0.0978300000000000
0.424750000000000	0.0969900000000000
0.464630000000000	0.0946300000000000
0.484550000000000	0.0931600000000000
0.504460000000000	0.0915100000000000
0.524350000000000	0.0897100000000000
0.544230000000000	0.0877500000000000
0.564100000000000	0.0856400000000000
0.583940000000000	0.0833600000000000
0.603770000000000	0.0809200000000000
0.623570000000000	0.0783200000000000
0.643360000000000	0.0755700000000000
0.663120000000000	0.0726700000000000
0.682860000000000	0.0696400000000000
0.702590000000000	0.0664700000000000
0.722290000000000	0.0631900000000000
0.741970000000000	0.0597700000000000
0.761630000000000	0.0562300000000000
0.781260000000000	0.0525500000000000
0.800870000000000	0.0487200000000000
0.820450000000000	0.0447600000000000
0.849760000000000	0.0385400000000000
0.869260000000000	0.0342100000000000
0.888730000000000	0.0297500000000000
0.908170000000000	0.0251400000000000
0.927570000000000	0.0204000000000000
0.946930000000000	0.0154900000000000
0.966250000000000	0.0104200000000000
1	0

I would like to be able to represent the curve created by these coordinates as a polynomial so that I can reproduce the curves to a decent degree of accuracy. This must be done as some downloaded airfoils do not contain equal numbers of upper and lower coordinates and I would like to be able to reproduce these airfoils with a user-defined number of x-coordinates ie. (0:0.01:1).

The two main concerns are as follows:

1. Input x coordinates are not evenly spaced so polyfit will skew the result as soon as I output it with evenly spaced x coordinates
2. First and final coordinates must not change so that any asymmetry/camber can be captured and polyfit fails to conform to this constraint.

Does anyone have any suggestions for a method here?

Answer 1

[John D'Errico]

xy = [0	0
    0.00369000000000000	0.0134500000000000
    0.0170200000000000	0.0280700000000000
    0.0336900000000000	0.0390000000000000
    0.0515300000000000	0.0479600000000000
    0.0699500000000000	0.0556700000000000
    0.107870000000000	0.0682100000000000
    0.127140000000000	0.0734700000000000
    0.146570000000000	0.0781200000000000
    0.166120000000000	0.0822000000000000
    0.185780000000000	0.0857400000000000
    0.205520000000000	0.0888000000000000
    0.225320000000000	0.0914200000000000
    0.245170000000000	0.0936200000000000
    0.265070000000000	0.0954000000000000
    0.284990000000000	0.0968000000000000
    0.304940000000000	0.0978000000000000
    0.324910000000000	0.0984500000000000
    0.344880000000000	0.0987500000000000
    0.364850000000000	0.0987300000000000
    0.384830000000000	0.0984200000000000
    0.404790000000000	0.0978300000000000
    0.424750000000000	0.0969900000000000
    0.464630000000000	0.0946300000000000
    0.484550000000000	0.0931600000000000
    0.504460000000000	0.0915100000000000
    0.524350000000000	0.0897100000000000
    0.544230000000000	0.0877500000000000
    0.564100000000000	0.0856400000000000
    0.583940000000000	0.0833600000000000
    0.603770000000000	0.0809200000000000
    0.623570000000000	0.0783200000000000
    0.643360000000000	0.0755700000000000
    0.663120000000000	0.0726700000000000
    0.682860000000000	0.0696400000000000
    0.702590000000000	0.0664700000000000
    0.722290000000000	0.0631900000000000
    0.741970000000000	0.0597700000000000
    0.761630000000000	0.0562300000000000
    0.781260000000000	0.0525500000000000
    0.800870000000000	0.0487200000000000
    0.820450000000000	0.0447600000000000
    0.849760000000000	0.0385400000000000
    0.869260000000000	0.0342100000000000
    0.888730000000000	0.0297500000000000
    0.908170000000000	0.0251400000000000
    0.927570000000000	0.0204000000000000
    0.946930000000000	0.0154900000000000
    0.966250000000000	0.0104200000000000
    1	0];
    
x = xy(:,1);
y = xy(:,2);
plot(x,y,'o')

So so often, people decide that a polynomial model is correct to use on their problem. Surely it must be! After all, people use polynomials all the time, so why would not it work here? Polynomials are so easy to use. And Taylor series are really just polynomials, and we learned to use them in school, so polynomials MUST be good, right?

Look at the plot I've drawn. Do you accept that the slope of that curve at x==0 will surely be essentially infinite? Or if you flip the axes around? When you do, the slope at that point will be zero, correct?

Again, at x=y=0, the curve when drawn as a function of x, it has a derivative singularity.

A feature of polynomial models is that they can NEVER represent a curve with infinite slope. Choose ANY polynomial, of ANY order. At no point on the finite real line will any polynomial EVER have an infinite derivative. This is a reason why for example, we cannot represent the function sqrt(x) as a Taylor series expanded around zero.

It does not matter how many terms you take in the expansion. There is NO polynomial model for this curve. Polyfit will never help.

Having said that, CAN you find a possibly valid model? For example, would a spline work? NO! A spline is built from cubic polynomial segments. Again, do cubic polynomials have an infinite slope? Do polynomials have singularities? NO. For example, suppose we fit a spline to that data? What is the derivative of the spline at x==0?

spl = spline(x,y);
fnval(fnder(spl),0)

ans = 4.5640

And that is a very finite number. The spline fit is not terrible. But you can do way better, with a much simpler to use model.

Ok, let me try again, is there a model you CAN use, that may work? Possibly, yes. The idea here would be to use a model that already has a singularity in a term, of the correct form.

For example. We have several fetures of your data that probably are important. We want the model to pass through the points (0,0) and (0,1). And we want it to have a derivative singularity at x==0, so an infinite derivative there. We might be thinking of a model in the form of a Pade approximant, for example. The nice thing about Pade approximations is they can have singularities in them. But the specific form we need is probably more specific yet. A Pade fit would not be exactly what we want.

Consider the model:

y(x) = sqrt(x)*(1-x)*P(x)

At x == 0, this curve will have a derivative singularity, AND at x==0, it will pass EXACTLY through zero. At x==1, the curve will pass EXACTLY thrrough zero again. It must at both endpoints, because sqrt(0)==0, and (1-1)==0.

All we need to do is now infer the function P(x). We might do that using least squares. I'll use the curve fitting toolbox here, although I could have done it using backslash too. (I chose the CFTB because the curve fitting TB allows me to do the fit simply, and it gives me nice plots.) And a 3rd or 4th degree polynomial for P(x) will be entirely adequate.

mdl = fittype('sqrt(x)*(1-x)*(a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4)','indep','x')

mdl = 
     General model:
     mdl(a0,a1,a2,a3,a4,x) = sqrt(x)*(1-x)*(a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4)

fittedmdl = fit(x,y,mdl)

Warning: Start point not provided, choosing random start point.

fittedmdl = 
     General model:
     fittedmdl(x) = sqrt(x)*(1-x)*(a0 + a1*x + a2*x^2 + a3*x^3 + a4*x^4)
     Coefficients (with 95% confidence bounds):
       a0 =      0.2117  (0.2103, 0.2131)
       a1 =      0.2461  (0.2277, 0.2645)
       a2 =     -0.4213  (-0.498, -0.3447)
       a3 =      0.2252  (0.1009, 0.3494)
       a4 =     0.04729  (-0.02074, 0.1153)

This is effectively a linear model in the unknown coefficients, so even though fit gets upset at me, I don't care. It will succeed without me providing starting values. And I'm lazy some days (ok most days.)

plot(fittedmdl,x,y)

So we have a very nice model, with only 4 coefficients needed. The model now has perfectly the properties we wanted. It has a derivative singularity at x==0. It passes EXACTLY through zero at x==0 and x==1. And it fits VERY nicely to the data. Do NOT forget that to achieve any accuracy, the coefficients of the cubic polynomial should be represented using more than 4 significant digits.

format long g
Pcoeffs = [fittedmdl.a4;fittedmdl.a3;fittedmdl.a2;fittedmdl.a1;fittedmdl.a0]

Pcoeffs = 5×1

0.0472869933538999 0.225152749466633 -0.421342589852467 0.246113708801069 0.211702199043612

So you could now use polyval to evaluate the model at any point too. Don't forget to multiply by those extra terms.

mdlfun = @(x) sqrt(x).*(1-x).*polyval(Pcoeffs,x);
mdlfun(0.5)

ans = 
        0.0921087676149181

I might expand the curve very carefully around zero.

xsmall = linspace(0,.11);
plot(xsmall,fittedmdl(xsmall),'r-',x(1:7),y(1:7),'bo')
xlim([0,.11])

As you can see, the curve now has accurately the shape you want it to have. Honestly, a cubic polynomial was probably entirely adequate too for P(X), as long as you carefully build the correct model.

Remember to look at your data. Always think about what properties it has. Can the model you have chosen possibly represent that curve? Polynomials cannot solve all problems. Although here, with just the correct nudging and a very careful case of model building, we found what is almost a polynomial, and it does perform very well indeed. (That sqrt(x) term makes it not technically a polynomial.)

[Alan Stevens]

@John D'Errico You're right, the spline behaviour near the zero point is bad!

[John D'Errico]

My suggestion on your answer fixes it though pretty nicely. It wants a singularity at that point. But we can put one into the curve. Which sort of makes the resulting spline a degree 3.5 spline. :)

Or, I suppose we could have carefully tweaked just the very first segment of the spline to have that singularity. But that would have forced the spline to be a hand made one.

[Alexander]

Thank you so much for your help, John. I've spent the best part of a few weeks tinkering with different solutions - iterating through Hicks-Henne-Bump and CST transform functions until Root Mean Square Error is tolerable but certainly didn't seem the most efficient, nor elegant solution so what appears, on the surface, to be a relatively simple problem.

I really appreciate the extra effort you have provided to help myself and other readers better understand your solution and you've provided plenty of context for future reading and learning.

I'll work on implementing this and let you know if it is a success!

[John D'Errico]

@Alexander

Thanks. That you spent some time trying to find a solution withut success just means you did not see the trick I employed to build a singularity directly into the model, and still have it easy to be fit. When you have built a gajillion models, that part gets easier. I'm happy it it helps you.

Answer 2

[Alan Stevens]

How about using a spline fit? After loading your data and calling the first column x and the second column y try

xfit = linspace(0,1,100);
yfit = spline(x,y,xfit);

plot(x,y,'.',xfit,yfit,'o-'),grid

[John D'Errico]

@Alan Stevens - see the answer I'm writing. Or if you want, just plot the shape of the curve to be fit. Then decide if a spline can represent that curve well. Because splines are still effectively polynomials, they also hate singularities.

[Alan Stevens]

That's exactly what I did!

xy = [0	0
    0.00369000000000000	0.0134500000000000
    0.0170200000000000	0.0280700000000000
    0.0336900000000000	0.0390000000000000
    0.0515300000000000	0.0479600000000000
    0.0699500000000000	0.0556700000000000
    0.107870000000000	0.0682100000000000
    0.127140000000000	0.0734700000000000
    0.146570000000000	0.0781200000000000
    0.166120000000000	0.0822000000000000
    0.185780000000000	0.0857400000000000
    0.205520000000000	0.0888000000000000
    0.225320000000000	0.0914200000000000
    0.245170000000000	0.0936200000000000
    0.265070000000000	0.0954000000000000
    0.284990000000000	0.0968000000000000
    0.304940000000000	0.0978000000000000
    0.324910000000000	0.0984500000000000
    0.344880000000000	0.0987500000000000
    0.364850000000000	0.0987300000000000
    0.384830000000000	0.0984200000000000
    0.404790000000000	0.0978300000000000
    0.424750000000000	0.0969900000000000
    0.464630000000000	0.0946300000000000
    0.484550000000000	0.0931600000000000
    0.504460000000000	0.0915100000000000
    0.524350000000000	0.0897100000000000
    0.544230000000000	0.0877500000000000
    0.564100000000000	0.0856400000000000
    0.583940000000000	0.0833600000000000
    0.603770000000000	0.0809200000000000
    0.623570000000000	0.0783200000000000
    0.643360000000000	0.0755700000000000
    0.663120000000000	0.0726700000000000
    0.682860000000000	0.0696400000000000
    0.702590000000000	0.0664700000000000
    0.722290000000000	0.0631900000000000
    0.741970000000000	0.0597700000000000
    0.761630000000000	0.0562300000000000
    0.781260000000000	0.0525500000000000
    0.800870000000000	0.0487200000000000
    0.820450000000000	0.0447600000000000
    0.849760000000000	0.0385400000000000
    0.869260000000000	0.0342100000000000
    0.888730000000000	0.0297500000000000
    0.908170000000000	0.0251400000000000
    0.927570000000000	0.0204000000000000
    0.946930000000000	0.0154900000000000
    0.966250000000000	0.0104200000000000
    1	0];
    
x = xy(:,1);
y = xy(:,2);
xfit = linspace(0,1,100);
yfit = spline(x,y,xfit);

plot(x,y,'.',xfit,yfit,'o-'),grid

Not perfect, but better than a polynomial fit.

[John D'Errico]

Not really. Again. LOOK VERY CAREFULLY NEAR ZERO. Do you see the strange bend there?

Check out my answer, as I have now finished it. It explains why a spline fails at zero. Yes, it does better than a normal polynomial. And I'll admit a spline model is decent there. And a spline does go exactly through the points, so that is a virtue.

[John D'Errico]

Ok. Let me suggest, that HAD you build a subtly different model, much of the form I built, a spline would have been perfectly adequate. For example, try this:

xy = [0	0
    0.00369000000000000	0.0134500000000000
    0.0170200000000000	0.0280700000000000
    0.0336900000000000	0.0390000000000000
    0.0515300000000000	0.0479600000000000
    0.0699500000000000	0.0556700000000000
    0.107870000000000	0.0682100000000000
    0.127140000000000	0.0734700000000000
    0.146570000000000	0.0781200000000000
    0.166120000000000	0.0822000000000000
    0.185780000000000	0.0857400000000000
    0.205520000000000	0.0888000000000000
    0.225320000000000	0.0914200000000000
    0.245170000000000	0.0936200000000000
    0.265070000000000	0.0954000000000000
    0.284990000000000	0.0968000000000000
    0.304940000000000	0.0978000000000000
    0.324910000000000	0.0984500000000000
    0.344880000000000	0.0987500000000000
    0.364850000000000	0.0987300000000000
    0.384830000000000	0.0984200000000000
    0.404790000000000	0.0978300000000000
    0.424750000000000	0.0969900000000000
    0.464630000000000	0.0946300000000000
    0.484550000000000	0.0931600000000000
    0.504460000000000	0.0915100000000000
    0.524350000000000	0.0897100000000000
    0.544230000000000	0.0877500000000000
    0.564100000000000	0.0856400000000000
    0.583940000000000	0.0833600000000000
    0.603770000000000	0.0809200000000000
    0.623570000000000	0.0783200000000000
    0.643360000000000	0.0755700000000000
    0.663120000000000	0.0726700000000000
    0.682860000000000	0.0696400000000000
    0.702590000000000	0.0664700000000000
    0.722290000000000	0.0631900000000000
    0.741970000000000	0.0597700000000000
    0.761630000000000	0.0562300000000000
    0.781260000000000	0.0525500000000000
    0.800870000000000	0.0487200000000000
    0.820450000000000	0.0447600000000000
    0.849760000000000	0.0385400000000000
    0.869260000000000	0.0342100000000000
    0.888730000000000	0.0297500000000000
    0.908170000000000	0.0251400000000000
    0.927570000000000	0.0204000000000000
    0.946930000000000	0.0154900000000000
    0.966250000000000	0.0104200000000000
    1	0];
x = xy(:,1);
y = xy(:,2);
spl = spline(x(2:end),y(2:end)./sqrt(x(2:end)));
fun = @(x) sqrt(x).*fnval(spl,x);
fplot(fun,[0,.1],'r-')
hold on
plot(x(1:7),y(1:7),'bo')

That would be a model that has the necessary behaviors at each end point, AND uses a spline. In fact, the resulting curve is actually better than mine in a sense, because it passes exactly through all of the data points, rather than being a least squares fit as I used. The important point is a spline does not perform well when you have a singularity in the curve.

Answer 3

[Bruno Luong]

You could use the free-knot spline

xy = [0	0
    0.00369000000000000	0.0134500000000000
    0.0170200000000000	0.0280700000000000
    0.0336900000000000	0.0390000000000000
    0.0515300000000000	0.0479600000000000
    0.0699500000000000	0.0556700000000000
    0.107870000000000	0.0682100000000000
    0.127140000000000	0.0734700000000000
    0.146570000000000	0.0781200000000000
    0.166120000000000	0.0822000000000000
    0.185780000000000	0.0857400000000000
    0.205520000000000	0.0888000000000000
    0.225320000000000	0.0914200000000000
    0.245170000000000	0.0936200000000000
    0.265070000000000	0.0954000000000000
    0.284990000000000	0.0968000000000000
    0.304940000000000	0.0978000000000000
    0.324910000000000	0.0984500000000000
    0.344880000000000	0.0987500000000000
    0.364850000000000	0.0987300000000000
    0.384830000000000	0.0984200000000000
    0.404790000000000	0.0978300000000000
    0.424750000000000	0.0969900000000000
    0.464630000000000	0.0946300000000000
    0.484550000000000	0.0931600000000000
    0.504460000000000	0.0915100000000000
    0.524350000000000	0.0897100000000000
    0.544230000000000	0.0877500000000000
    0.564100000000000	0.0856400000000000
    0.583940000000000	0.0833600000000000
    0.603770000000000	0.0809200000000000
    0.623570000000000	0.0783200000000000
    0.643360000000000	0.0755700000000000
    0.663120000000000	0.0726700000000000
    0.682860000000000	0.0696400000000000
    0.702590000000000	0.0664700000000000
    0.722290000000000	0.0631900000000000
    0.741970000000000	0.0597700000000000
    0.761630000000000	0.0562300000000000
    0.781260000000000	0.0525500000000000
    0.800870000000000	0.0487200000000000
    0.820450000000000	0.0447600000000000
    0.849760000000000	0.0385400000000000
    0.869260000000000	0.0342100000000000
    0.888730000000000	0.0297500000000000
    0.908170000000000	0.0251400000000000
    0.927570000000000	0.0204000000000000
    0.946930000000000	0.0154900000000000
    0.966250000000000	0.0104200000000000
    1	0];
x = xy(:,1);
y = xy(:,2);

pp=BSFK(x,y); % https://fr.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation

foilfun = @(x) ppval(pp,x);
xi = linspace(0,1,501);

close all
plot(x,y,'or',xi,foilfun(xi),'b-')