# How to represent [airfoil] coordinates as a polynomial

27 views (last 30 days)

Show older comments

Alexander
on 12 Jan 2023

Commented: John D'Errico
on 12 Jan 2023

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:

- Input x coordinates are not evenly spaced so polyfit will skew the result as soon as I output it with evenly spaced x coordinates
- 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?

##### 0 Comments

### Accepted Answer

John D'Errico
on 12 Jan 2023

Edited: John D'Errico
on 12 Jan 2023

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)

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')

fittedmdl = fit(x,y,mdl)

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]

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)

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.)

##### 4 Comments

John D'Errico
on 12 Jan 2023

### More Answers (2)

Alan Stevens
on 12 Jan 2023

Edited: Alan Stevens
on 12 Jan 2023

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

##### 4 Comments

John D'Errico
on 12 Jan 2023

Edited: John D'Errico
on 12 Jan 2023

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.

Bruno Luong
on 12 Jan 2023

Edited: Bruno Luong
on 12 Jan 2023

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-')

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!