2nd order relationship between density and altitude using curve fitting.

Please someone help. I'm trying this for many hours but couldn't get desired result.Actually I'm trying to obtain the 2nd order relationship between density vs altitute.I used different methods like polyfit,fit,fitlm,etc but didn't succeed. I need to calculate the density at 7000m that is 7km and answer(density) should be around 0.6kg/m^3 but my answer is coming very much around 94000kg/m^3.There are many solutions of this problem are available but not with MATLAB(They are EES,etc).I could not understand how did they achieved the solution. Here is my work

 clc;clear;close all;
%Data
z = [6377 6378 6379 6380 6381 6382 6383 6385 6387 6392 6397 6402]; %km
rho = [1.225 1.112 1.007 0.9093 0.8194 0.7364 0.6601 0.5258 0.4135 0.1948 0.08891 0.04008];  %kg/m^3
%2nd order Relationship between rho and z using polyfit
coeff = polyfit(z,rho,2)
%evaluate the density at z using this relation
rho_cal = polyval(coeff,z);
%plotting
plot(z,rho,'*',z,rho_cal)
xlabel('z [km]'); ylabel('\rho [kg/m^3]')
title('Atmospheric density with elevation')
legend('Data','Relation')
%density at 7000m altitude that is 7km
 rho = polyval(coeff,7)

Also they achieved the polynomial fitted curve solution rho=0.00223747z^2-0.10167z+1.20252 please help me.I couldn't figure out what I'm missing.

Answer 1

[John D'Errico]

z = [6377 6378 6379 6380 6381 6382 6383 6385 6387 6392 6397 6402]; %km
rho = [1.225 1.112 1.007 0.9093 0.8194 0.7364 0.6601 0.5258 0.4135 0.1948 0.08891 0.04008];  %kg/m^3

That is the data. z is in km, and it goes no at all as high as 7000km, but you apparently want to extrapolate out that far. First, PLOT YOUR DATA. ALWAYS PLOT EVERYTHING.

plot(z,rho,'o')

Ok. Now you want to extrapolate that out as far as 7000. But you used a quadratic polynomial. SERIOUSLY? Really, be serious, and think about what shape a quadratic polynomial would have. First, I'll fit a quadratic to your data. Since you used polyfit, I'll do that.

P2 = polyfit(z,rho,2)

P2 = 1×3

1.0e+04 * 0.0000 -0.0029 9.1639

rho7000 = polyval(P2,7000)

rho7000 = 806.2896

So the quadratic polynomial prediction out that far is over 800. But again, think about the fact that a quadratic polynomial has the shape of a parabola. It will certainly turn around, and start to shoot upwards. Even if I extrapolate out only as far as 6500, look at the shape you see.

zpred = linspace(6375,6500);
rhopred = polyval(P2,zpred);

plot(z,rho,'o',zpred,rhopred,'r-')

By the time it gets out as far as 7000, what do you expect? You have data over a terribly narrow rangee, and you want to extrapolate out where? Do I need to post my favorite quote by Mark Twain and the dangers of extrapolation?

Mark Twain, Life on the Mississippi (1884):

“In the space of one hundred and seventy six years the Lower Mississippi has shortened itself two hundred and forty-two miles. That is an average of a trifle over a mile and a third per year. Therefore, any calm person, who is not blind or idiotic, can see that in the Old Oölitic Silurian Period, just a million years ago next November, the Lower Mississippi was upwards of one million three hundred thousand miles long, and stuck out over the Gulf of Mexico like a fishing-pole. And by the same token any person can see that seven hundred and forty-two years from now the Lower Mississippi will be only a mile and three-quarters long, and Cairo [Illinois] and New Orleans will have joined their streets together and be plodding comfortably along under a single mayor and a mutual board of aldermen. There is something fascinating about science. One gets such wholesale returns of conjecture out of such a trifling investment of fact.”

A good rule is never extrapolate far beyond the support of your data. How far is far? Sigh, that kind of depends on how noisy is your data, how nonlinear is it, and what model I was using to extrapolate. Does the model make physical sense? That is, some models are used that are based physical reasonaing about the data and where the relationship comes from. These models might be trusted more, than say a polynomial model. A polynomial model has no reason to be correct, or to be trusted beyond the region of support.

Now maybe you were told to use such a model by your teacher. How can we know?

Just looking at your data, it appears to be something that is approaching an asymptote, perhaps eventually flattening out. So you might use a model that has that property. Such a model may have been chosen for no good reason of course, but it MAY be a better choice than a quadratic. For example, better might be to use a simple negative exponential model. Something of this form:

rho = a*exp(-b*z)

If that is the case, then a semilogy plot will be a straight line.

semilogy(z,rho,'-o')

So not terrible, but not terribly great either. But at least it will predict something rational if you extrapolate it. And that model is easy to generate from polyfit. Just log z, then use a linear model.

P1 = polyfit(z,log(rho),1)

P1 = 1×2

-0.1365 870.9080

exp(polyval(P1,7000))

ans = 1.6159e-37

so it predicts a value that is essentially zero if we go out that far. And that may be unreasonably low, but at least it is far more meaningful that 806. Actually, you might have decided to use a quadratic polynomial for log(rho). At least that polynomial will not turn around and start predicting completely meaningless results, because this data has negative curvature, even in log(rho).

P2 = polyfit(z,log(rho),2)

P2 = 1×3

1.0e+04 * -0.0000 0.0022 -6.8786

exp(polyval(P2,7000))

ans = 3.9890e-314

As expected, that predicts a value even far closer to zero when you go that far out. Still better than 806 though.

[ajeet sahu]

Firstly thank you very much john for answering it in a detailed way.Actually this is standard textbook problem.

I'd checked the solution of this Problem online on different platforms(chegg,quizlet,etc) even the official solution manual but everywhere the solution was same.Means they fitted the parabola using EES,etc and their solution was rho=0.00223747z^2-0.10167z+1.20252 and then they just put z=7(7000m means 7km) to get the density of 0.6kg/m^3.I was wondering how did they achieved this by using a parabola but why I couldn't.Now I've plotted their fitted parabola and I'm sure that solution was wrong. I wasted many hours to prove(get) something that doesn't exist. Thank you again john...

[John D'Errico]

Let me see if I understand.

They gave you data that has a support of between [6377,6402], they ask you to predict a result at z = 7000. This is perhaps the most dumba$$ problem writer I have ever seen, if this is what they expect you to predict. Doing that with a quadratic polynomial is flat out insanity.

As I said, the only reason I could imagine doing that is if the goal was to teach a student why to NOT do that. So I might hope they go into a detailed explanation as to why it is such a terribly bad idea.

Anyway, the model you show does not even remotely come close to predict the data, regardless if I scale z by dividing by 1000 or not. Try plotting that polynomial over the region of your data. Soes it predict rho, even if you are not extrapolating? It seems to not come even close, regardless if I scale z by dividing by 1000.

rho_z = @(z) 0.00223747*z.^2-0.10167*z+1.20252;

rho_z(7000)

ans = 1.0893e+05

Clearly that is meaningless, so I presume this model uses z in km.

rho_z(7)

ans = 0.6005

which does evaluate to 0.6005. So one presumes they are scaling by 1000.

z = [6377 6378 6379 6380 6381 6382 6383 6385 6387 6392 6397 6402]; %km
rho = [1.225 1.112 1.007 0.9093 0.8194 0.7364 0.6601 0.5258 0.4135 0.1948 0.08891 0.04008];  %kg/m^3

zpred = linspace(6.370,7,1000);
plot(z/1000,rho,'bo',7,rho_z(7),'gx',zpred,rho_z(zpred),'-r')

So the red line is the curve they suggest is the correct model to predict that data, and to extrapolate that relationship? SERIOUSLY?

Even over the support of the data, we see this:

zpred = linspace(6.370,6.402);
plot(z/1000,rho,'bo',zpred,rho_z(zpred),'-r')

That data is actually pretty smooth. It shows a significant amount of signal, and almost no noise. So I would be quite happy modeling it with some simple function on that support. Even a quadraitc polynomial is fine over that short of a domain. Just don't extrapolate it, at least not far.

I'd look for a better textbook. I'm not one for burning books, but if I was cold in the winter, that book would be first on my list to use for kindling.