How to represent [airfoil] coordinates as a polynomial

33 ビュー (過去 30 日間)
Alexander
Alexander 2023 年 1 月 12 日
コメント済み: John D'Errico 2023 年 1 月 12 日
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?

採用された回答

John D'Errico
John D'Errico 2023 年 1 月 12 日
編集済み: John D'Errico 2023 年 1 月 12 日
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.)
  4 件のコメント
Alexander
Alexander 2023 年 1 月 12 日
移動済み: Bruno Luong 2023 年 1 月 12 日
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
John D'Errico 2023 年 1 月 12 日
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.

サインインしてコメントする。

その他の回答 (2 件)

Alan Stevens
Alan Stevens 2023 年 1 月 12 日
編集済み: Alan Stevens 2023 年 1 月 12 日
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 件のコメント
John D'Errico
John D'Errico 2023 年 1 月 12 日
編集済み: John D'Errico 2023 年 1 月 12 日
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
John D'Errico 2023 年 1 月 12 日
編集済み: John D'Errico 2023 年 1 月 12 日
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
Bruno Luong 2023 年 1 月 12 日
編集済み: Bruno Luong 2023 年 1 月 12 日
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-')

カテゴリ

Help Center および File ExchangeSpline Postprocessing についてさらに検索

製品


リリース

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by