gauss elimination or inverse

Hi all,

I have a plynomial equation similar to the one iven below

Equation : ax^(n-1)+bx^(n-2)+.............+c=y

The equation must be satisfied for a set of x (say x1,x2,x3.....etc) and correspondig y (say y1,y2,y3......etc)

I need a suggestion to find the coefficient for the equation.(a,b,.....,c)

Should i use Gauss elimination method or simply do a inverse for the matrix [x] and find the solution.

Highly appreciate the responses.

Answer 1

[John D'Errico]

People have said repeatedly that you want to use polyfit. This is the correct answer, IF you need to fit a polynomial.

At the same time, you should not use Gaussian elimination, or use a matrix inverse. Those tools will have you writing the code, and sadly, making many novice mistakes in the computations that are already dealt with in polyfit.

But even with polyfit, there are still massively serious issues when you are trying to fit a high order polynomial. These issues can make it almost impossible to do a high order polynomial fit in double precision arithmetic. The problem is the resulting linear system of equations will become nearly singular in double precision arithmetic. (It would be far worse had you tried to use Gaussian elimination or use a matrix inverse. Seriously so.)

The big problem is what domain the variable x lives in. What is the min and max for x? You don't show your data. How many datapoints do you have? How much noise is present in the data? Best, is if you provide your data. Use a .mat file, attached to a comment.

In order to have even a chance in hell of solving this problem, you need to use the centering and scaling options in polyfit. That will improve the condition of the solution. But even then, unless we know enough about your data, it is difficult to help you more, to know if you can actually do this fit. There is potentially a chance in the end that you can do nothing at all to solve the problem, but I cannot know that until I see your data.

Finally, I would ask if you really do need a polynomial to fit the data? High order polynomial fits are usually a terrribly bad idea, when far better tools, like spline models, are available.

[deep]

Attached are the data points for the curve.

Any suggestion on the same is highly appreciated.

Hoping for something to work out.

Thank you for the response.

[John D'Errico]

Ok. This is one of the things I was worrried about. There are actually several problems, although one is a complete show stopper. You data has several variables in it. One of them is z. What is z? But you are asking about modeling x and y.

[x,y]
ans =
            1            1
         1.01         2.13
       1.0078        3.261
       1.1205        4.385
       1.3478        5.492
        2.182        7.586
        2.817         8.52
         3.58        9.353
        4.429       10.099
        6.265        11.42
        7.225        12.01
        8.203        12.58
        9.188        13.13
        10.18        13.68
        12.19        14.72
        14.24        15.66
        16.36        16.44
        18.56        16.95
        20.82         17.1
plot(x,y,'o-')
grid on

Do you see a major problem here? What polynomials can you think of that have an infinite slope? Thus a singularity? You CANNOT fit this curve with a direct polynomial model. No way. No possibility. Not ever. Polynomials do not have a singularity.

It depends on what you want to do. Are you looking for a function that you can write down? Do you just want to create a smoooth curve that goes through the points for a nice plot? Whatis your real need?

Can you do anything at all? Hmm. well,first, lets go back up and look at your data.

1            1
1.01         2.13
1.0078        3.261
1.1205        4.385
...

So the third data point is not in any useful order. It looks like x(3) is wrong. The x values are not in order, yet y(3) is right in the middel of y(2) and y(4).So your data has some garbage in it. Probably transcription error is my guess. I would wonder if x(3) should have been 1.078, not 1.0078.

Yet, can we do ANYTHING? A nice trick is to see if this curve can be transformed for a fit. The idea is to find a model to y(x) in the end.

plot(nthroot(x([1 2 4 5 6 7])-1,3),y([1 2 4 5 6 7]),'-o')
grid on
xlabel 'Cube root of x'
ylabel 'y'

So we see that a cube root transformation on (x-1) results in a nice linear curve down near zero. Essentailly, I was able to remove that essential singularity. Note that I just dropped out that third garbage point here. See that sqrt and 4th roots have the wrong shape down near zero.

plot(nthroot(x - 1,3),y,'-o')
grid on
xlabel 'Cube root of x'
ylabel 'y'

As I said, x(3) is garbage. But now, if we try to model y as a function of the cube root of (x-1), we can at least have some tiny chance of success.

ind = [1 2 4:19];
format long g
P8 = polyfit(nthroot(x(ind)-1,3),y(ind),8);
P8'
ans =
         0.665189276257975
         -7.95919647163279
          38.7832233452933
         -99.3344261623917
          143.598072175124
         -115.901599470974
          46.9978270959549
        -0.631331340664469
         0.996380665885126

an 8th degree polynomial fit is actually very near the highest I would ever recommend. But this seems to have done reasonably well.

As you should expect, things go completely to hell for x even slightly below x==1 or above x==20.82. Remember, this model is a model in the form of:

y = f(nthroot(x-1,3))

So you have a polynomial in powers of (x-1)^(1/3). Even at that, splines would still be my ultimate preference. And again, you could never in a million years have fit a true polynomial in pure integer powers of x itself.

Finally, looking at the data, I wonder if you might know the original (untransformed) curve should go EXACTLY through the point (1,1)? That can be achieved easiy enough.

[deep]

Hi John,

I highly appreciate your elaborate explanation on the doubt. It clarifies a lot of information and also insight to what i had have missed out.

What is z? :

• here z is x, file was not renamed as z

Are you looking for a function that you can write down? Do you just want to create a smoooth curve that goes through the points for a nice plot? Whatis your real need?

• Yes, I need to generate smooth curve through the given points for various stations of a hull.These are offset points for the same.My real need is to generate curves that pass through set of such numerical data (offsets).

• If you see plot [y,x] you will see a section of the hull form.

• I am aware of the singularity problem arising in the given data, although an expert advice on the same was required as i have limited knowledge on matlab capabilities.

Finally, looking at the data, I wonder if you might know the original (untransformed) curve should go EXACTLY through the point (1,1)? That can be achieved easiy enough.

• This is the original curve. Achieving the shape of the curve is important here, this is just a sample for one such curve and i have many such curves to be fitted. The base point for each must be same may be say(0,0) or (1,1) or any other point. As all these curves lie longitudinally (here the z axis) on the same base point (y,x) (or points if the shape changes i.e. the base flattens)

[John D'Errico]

From your comments, you have no need for a polynomial form at all. All you seem to be looking to achieve is a fit that passes smoothly through the points. So there is not even a remote need to use polynomial models, or even the contrivance I came up with to model

Not surprisingly, this is what splines do very well. In fact, the word spline (lath is also used in this respect) can be found in ship building, a tool used to draw the smooth curves needed when lofting hull shapes.

https://en.wikipedia.org/wiki/Flat_spline

You will still have a problem trying to interpolate this curve using the standard spline tools in MATLAB, because they are not designed to fit curves with derivative singularities in them. (Which is what you have.) Instead, you want to use a tool like my interparc. It can be freely downloaded from the file exchange, from this link:

https://www.mathworks.com/matlabcentral/fileexchange/34874-interparc

Use of interparc is trivial once you download it. I don't even need to exclude the third point from this interpolation, even though I know the value of x was an error. That is because interparc has no problems with that tiny deviation. It still draws a smooth curve through the points.

XY = interparc(100,x,y);
plot(XY(:,1),XY(:,2),'-',x,y,'o')

So 100 points smoothly interpolated along that curve. They are chosen to be equidistant from each other along that curve. And the curve will pass exactly through those points.

So the answer to your problem is to NOT use any polynomial at all. Instead, just use interparc.

[deep]

Thank you for your elaborate explanation.

Answer 2

[Stephan]

Hi,

use polyfit to do this.

Best regards

Stephan

[deep]

Thank you,

I applied the polyfit method to my set of data and it does not result in the correct results as the order of my equations are expected to be higher than 10.

for eg:

y1=ax1^15+bx1^14+......+C1

y2=ax1^15+bx1^14+......+C2

.

.

.

.

.

Y16=ax16^15+bx16^14+......C16

How to work out this to get the coefficients correct?

Thank you in advance

[Bruno Luong]

How to work out this to get the coefficients correct?

Give up such high order fitting. Change the method.

[deep]

This is my first look on the idea.

I will try to look for suitable alternative method.

Thank you for your advice.

Answer 3

[Bruno Luong]

P = polyfit(x(:),y(:),n-1);
a = P(1);
b = P(2);
...
c = P(n);

[deep]

Will this work for orders as high as 25 i.e. n=25

[Bruno Luong]

In general fitting polynomial with order >= 10 is not reliable.

[deep]

I applied the polyfit method to my set of data and it does not result in the correct results as the order of my equations are expected to be higher than 10.

for eg:

y1=ax1^15+bx1^14+......+C1

y2=ax1^15+bx1^14+......+C2

.

.

.

.

.

Y16=ax16^15+bx16^14+......C16

How to work out this to get the coefficients correct?

Thank you in advance

Answer 4

[Walter Roberson]

polyfit(x(:), y(:), n - 1)

[deep]

I applied the polyfit method to my set of data and it does not result in the correct results as the order of my equations are expected to be higher than 10.

for eg:

y1=ax1^15+bx1^14+......+C1

y2=ax1^15+bx1^14+......+C2

.

.

.

.

.

Y16=ax16^15+bx16^14+......C16

How to work out this to get the coefficients correct?

Thank you in advance

[Torsten]

What do you mean by "it does not result in the correct results" ?

[deep]

The coeffients do not match the values obtained from manual calculations.

I.e. if i get the values of a,b,c........ from polyfit

I use them to calculate the y1 using the x1 value and the coefficients

y1=ax1^15+bx1^14+......+C1

So here theres a mismatch observed.

[Torsten]

If you insist on fitting your data with a polynomial of degree 15 (which is numerical nonsense in my opinion), use the centering and scaling option of "polyfit" as described here:

https://de.mathworks.com/help/matlab/ref/polyfit.html

[deep]

This is my first look on the idea.

I will try to look for suitable alternative method.

Thank you for your advice.