Thread Subject:

Problem

My problem is following:

I have velocity u and depth z data matrix (1x100). And have to find answer to the du/dz.

I know the answer when I have u like a equation - u=2*z+21. But how i manage this problem when I have numerical input data.

Thanks

in Matlab

Beaver <marden.nomm@gmail.com> wrote in message <1892926589.71940.1273005012321.JavaMail.root@gallium.mathforum.org>...
> in Matlab

that's exactly right...

us

Beaver wrote:
> My problem is following:
>
> I have velocity u and depth z data matrix (1x100). And have to find answer to the du/dz.
>
> I know the answer when I have u like a equation - u=2*z+21. But how i manage this problem when I have numerical input data.

It is impossible to do when you have a finite amount of sampled input data; it
would even be impossible to do if you had an continuous (infinite) amount of
data that had been truncated to a finite precision. You cannot even always do
it if you have a continuous infinite amount of infinitely precise data.

In order to arrive at du/dz (which is a quantity defined for The Calculus Of
The Infinitesimals, and thus inherently requires infinite representation), you
need a symbolic function which is continuous -- though of course if the
symbolic function was finitely piecewise continuous, you could come up with a
finitely piecewise differential that was valid except at the transition points.

If you have a finite sampling of data, there are _always_ an infinite number
of different functions that might be the "real" function, with all of those
infinite number of functions exactly matching the known data at the known
locations.

In such a situation, if you have a model of what form the function takes, and
the function has a number of parameters with unknown values with the number of
parameters not exceeding the number of data samples you have, then you _might_
be able to calculate the unknown parameters -- but not usually, not unless you
have at most 2 data points and the parameters happen to be exact multiples of
powers of 2... You might be able to _approximate_ the parameters -- that is
somewhat more likely... if the model happens to take a fairly restricted form
such as being a pure polynomial...

Of course, your success is going to depend upon having chosen the correct
model to start with -- and most people don't. For example, there is a tendency
for people who have K data points to model their data as a polynomial of
degree K. They get out numeric solutions that might match the data as closely
as one might reasonably hope for... and yet which will very likely be the
wrong polynomial. For example, if I give you 100 points that sample what is
really a cubic polynomial, then constructing a degree-100 polynomial is almost
certainly _not_ going to result in the coefficients for degrees 4 to 100 to be
found to be precisely zero... they might be found to be quite small, but as
far as the program is concerned, there is no reason to know that the
calculated values are "noise" or genuine coefficients.


Science is a cruel task-master. It is not uncommon in science to develop a
"working model" of what is happening, and to dig into it for *years* and come
up with theories... and then a decade or a century down the road, for someone
to find new evidence or a new way of thinking about the science that shows
that a notably different mathematical form is better.

For example, the data that you have for river depth -- do you also have
readings for the solar flux that was pushing on the surface of the water, and
thus altering the depth? Sure, the solar flux might only have an effect of
(say) one part in 60,000 -- but when you are calculating a 100 degree
polynomial, a measurement error of 1 part in 60000 is enough to cause the
second coefficient (for degree 99) to change by a factor of 600, and the last
coefficient (the constant term) to change by a factor of 1 in 6E477 . Even if
you manage to avoid numeric overflow in the calculations, that would be about
like changing the constant term from 6E+235 to 1E-235 . Is your data really
measured that accurately? If not, then clearly any answers you got out from
fitting with a 100 degree polynomial would be physical nonsense.

On May 5, 8:29 am, Beaver <marden.n...@gmail.com> wrote:
> My problem is following:
>
> I have velocity u and depth z data matrix (1x100). And have to find answer to the du/dz.
>
> I know the answer when I have u like a equation - u=2*z+21. But how i manage this problem when I have numerical input data.
>
> Thanks

This is the absolute core of calculus.
Do you remember in High School when you learned about derivatives?
What is a derivative exactly?
Clue: it is a limit.

In numerical analysis, we cannot take limits, so we must approximate.
That's enough clues.
Come back to us with a Matlab question.

BTW, a variable that is (1x100) is not a matrix. It is a row vector.

Beaver <marden.nomm@gmail.com> wrote in message <1606833116.71939.1273004986245.JavaMail.root@gallium.mathforum.org>...
> My problem is following:
>
> I have velocity u and depth z data matrix (1x100). And have to find answer to the du/dz.
>
> I know the answer when I have u like a equation - u=2*z+21. But how i manage this problem when I have numerical input data.
>
> Thanks
- - - - - - - -
  As you can gather from the previous remarks, you can only obtain an approximate derivative from discrete data. The closeness of the approximation will depend on the fineness and accuracy of the data.

  If you had only two points (z1,u1) and (z2,u2), your best guess would be the straight line slope between the points

 (u2-u1)/(z2-z1)

and from this you can extend this idea using matlab's 'diff' function to

 diff(u)./diff(z)

for longer vectors. Unfortunately this rather crude method gives one less derivative approximation than the number of values in the vectors, since each one represents a derivative approximation for a point midway between data points, and also it gives only a linear type of approximation.

  I will now give you a formula that produces as many derivative approximations as in the given vectors with each one centered over a corresponding point in those vectors, and which in addition amounts to a second order approximation. That is, if the one vector's values are a second order polynomial function of values in the other vector, these derivatives would be exact (save for round off error, of course.) I will give it in terms of variables x and y rather than your z and u. If x and y are the given vectors (of the same length,) do this:

 x1 = x([3,1:end-1]); x2 = x([2:end,end-2]);
 y1 = y([3,1:end-1]); y2 = y([2:end,end-2]);
 dydx = ((y2-y).*(x-x1).^2+(y-y1).*(x2-x).^2)./ ...
        ((x2-x).*(x-x1).*(x2-x1));

The 'dydx' vector will be an approximation to the derivative dy/dx. One requirement for using this formula is that the vectors must possess at least three points.

  There are of course higher order formulas with correspondingly more complicated formulas. I trust this one will be adequate for your purposes.

Roger Stafford

Thanks to all. I find good advices.

Is this formula have some notes how this works?

Beaver <marden.nomm@gmail.com> wrote in message <1317322228.78895.1273096933766.JavaMail.root@gallium.mathforum.org>...
> Thanks to all. I find good advices.
>
> Is this formula have some notes how this works?
- - - - - - -
  Prepare yourself for a long-winded explanation. If we assume that the x values are increasing, then, except for the endpoints, the quantities xa=x1(i), xb=x(i), and xc=x2(i) at each i-th point satisfy xa < xb < xc, and we are endeavoring to approximate the derivative, dy/dx, at the intermediate xb value.

  Defining ya=y1(i), yb=y(i), and yc=y2(i), we do this by determining the quadratic function,

 Y = A*(X-xb)^2+B*(X-xb)+C,

which goes through the three points (xa,ya), (xb,yb), and (xc,yc). Clearly C must be yb. Finding the other two coefficients, A and B, for this quadratic involves solving two equations in the two coefficient unknowns. The derivative of that quadratic function is

 dY/dX = 2*A*(X-xb)+B,

and when X is set to xb, this is simply B. After the equations are solved, this derivative value, B, turns out to be

 B = ((yc-yb)*(xb-xa)^2+(yb-ya)*(xc-xb)^2)/((xc-xb)*(xb-xa)*(xc-xa))

which agrees with the formula given earlier in this thread.

  Actually the assumption xa < xb < xc is not needed in the above formula, and that fact is used in the case of the two endpoints. By the addition of the third point in [3,1:end-1] to the left end of x1 and y1, and the "end-2" point in [2:end,end-2] to the right end of x2 and y2 we have

 xb < xc < xa

in the first case and

 xc < xa < xb

in the second case at these endpoints. Thus in these cases the derivative is being evaluated at the side of the interval containing the three points rather than at the point in its interior. Nevertheless the same formula remains valid for finding the derivative of the quadratic even though now being evaluated at a side point. That is why the endpoints do not have to be treated as a special case in this method.

  Another way to envision this formula is to write it in this equivalent form:

 dydx = ((xb-xa)/(xc-xa)) * ((yc-yb)/(xc-xb)) + ...
        ((xc-xb)/(xc-xa)) * ((yb-ya)/(xb-xa))

The quotients on the right are the two straight-line derivatives between points b and c and between a and b. The two quotient factors on the left involving only x-coordinates can be considered weights applied to these two straight-line derivatives, with their sum always equal to one. In other words the formula in question is nothing more than a certain weighted average between the two straight-line derivatives on either side of the point (xb,yb).

  I hope this explanation will indicate the reasoning behind the formula I gave you.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hrt7e0$drs$1@fred.mathworks.com>...
> Prepare yourself for a long-winded explanation. .......
- - - - - - - - -
  Please add this onto the previous discussion I gave.

  If y = f(x) is some continuous function with its first three derivatives also continuous over the interval xa < xb < xc, then by an extension of the famous Rolle's theorem of elementary calculus, there is a value xi (Greek letter) such that the difference between the expression above for an approximation of the derivative at xb and the actual derivative, f'(xb) is equal to

  1/6*f'''(xi)*(xc-xb)*(xb-xa)

where xi lies in xa < xi < xc and f'''(xi) is the third derivative value of f(x) there. A remarkable fact!

  What this signifies is that if the function whose derivative you are attempting to approximate with the formula I gave you has a third derivative which remains very small - in other words is reasonably smooth - then you are guaranteed a close approximation.

Roger Stafford

On May 6, 8:15 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
> "Roger Stafford" <ellieandrogerxy...@mindspring.com.invalid> wrote in message <hrt7e0$dr...@fred.mathworks.com>...
> >   Prepare yourself for a long-winded explanation. .......
>
> - - - - - - - - -
>   Please add this onto the previous discussion I gave.
>
>   If y = f(x) is some continuous function with its first three derivatives also continuous over the interval xa < xb < xc, then by an extension of the famous Rolle's theorem of elementary calculus, there is a value xi (Greek letter) such that the difference between the expression above for an approximation of the derivative at xb and the actual derivative, f'(xb) is equal to
>
>   1/6*f'''(xi)*(xc-xb)*(xb-xa)
>
> where xi lies in xa < xi < xc and f'''(xi) is the third derivative value of f(x) there.  A remarkable fact!
>
>   What this signifies is that if the function whose derivative you are attempting to approximate with the formula I gave you has a third derivative which remains very small - in other words is reasonably smooth - then you are guaranteed a close approximation.
>
> Roger Stafford

But Roger, doesn't your fancy algorithm reduce to central finite
differences if the x are equispaced:
dydx(2:end-1)=(y(3:end) - y(1:end-2))/(2*dx);
without the pretty stuff at the ends of course.

TideMan <mulgor@gmail.com> wrote in message <c29d4c47-5dad-44ea-879c-f3e5c682148d@o8g2000yqo.googlegroups.com>...
> But Roger, doesn't your fancy algorithm reduce to central finite
> differences if the x are equispaced:
> dydx(2:end-1)=(y(3:end) - y(1:end-2))/(2*dx);
> without the pretty stuff at the ends of course.

  Yes, of course, for equally-spaced x values. In that case one can obtain the same answers from the single-dimensional version of the 'gradient' function except for the endpoints, which are not properly done in my opinion.

  For variable-spaced x's the 'gradient' function does not use weighted averages, so it is not actually a second order approximation in that circumstance. That is, it does not get exact derivatives for quadratic functions if the x intervals are variable - at least not in my version.

Roger Stafford

Thanks, its works perfectly

By the way is it formula has a name or something. I have to mark this to my work.

Marden

Beaver <marden.nomm@gmail.com> wrote in message <182806353.128131.1273733858816.JavaMail.root@gallium.mathforum.org>...
> Thanks, its works perfectly
>
> By the way is it formula has a name or something. I have to mark this to my work.
>

Savitzky-Golay filter.

http://www.mathworks.com/access/helpdesk/help/toolbox/signal/sgolay.html

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hsg93c$kkk$1@fred.mathworks.com>...
> Beaver <marden.nomm@gmail.com> wrote in message <182806353.128131.1273733858816.JavaMail.root@gallium.mathforum.org>...
> > Thanks, its works perfectly
> > By the way is it formula has a name or something. I have to mark this to my work.
> > Marden
>
> Savitzky-Golay filter.
> http://www.mathworks.com/access/helpdesk/help/toolbox/signal/sgolay.html
> Bruno

To Bruno:

  It would only be considered a Savitzky-Golay filter if the x-coordinates of the data were equally-spaced. However, the whole point of my suggestion was to handle the case of unequally-spaced points.

 http://en.wikipedia.org/wiki/Numerical_smoothing_and_differentiation
 http://en.wikipedia.org/wiki/Savitzky–;Golay_smoothing_filter

To Marden:

  This formula can be derived directly from Lagrange's interpolation formula for three points by evaluating the derivative at the central point.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hsgls4$7gh$1@fred.mathworks.com>...
>
>
> This formula can be derived directly from Lagrange's interpolation formula for three points by evaluating the derivative at the central point.

Roger, do you have any idea about the robustness of such method when the point abscissa get closer? It's surely not good, but is it more or less robust than the diff command (order=1)? Or more generally what happens when the polynomial order increases?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hsgp6v$d9d$1@fred.mathworks.com>...
> Roger, do you have any idea about the robustness of such method when the point abscissa get closer? It's surely not good, but is it more or less robust than the diff command (order=1)? Or more generally what happens when the polynomial order increases?
> Bruno

  I don't claim to be an expert in the subject of noise and filters, so I'll just do a little "hand waving".

  With data that is comparatively free of noise, the question of low order versus high order polynomial approximations for finding derivatives is dependent on the relative sizes of the higher order versus lower order derivatives inherent in the data. Where higher order derivatives are comparatively small, then high order polynomial approximation is called for, but low order polynomials are better when the data exhibits irregular behavior in its higher derivatives - in others words when it is not so "smooth".

  In the absence of noise, close spacing is of course all to the good, the closer the better, especially for the higher order polynomials. The mean value theorem's error term contains not only a high order derivative but is multiplied by a corresponding high power of the spacing interval which one would like to see as small as possible.

  However, when noise enters the picture in a serious way, everything changes. Close spacing now becomes undesirable and high order polynomials tend to become worse than low order. But even low order derivative methods suffer from the noise and increasingly require larger spacing as the noise increases. Therefore anyone trying to find derivatives from discrete data needs to take careful account of just how noisy that data is likely to be in choosing the best method.

  On the question of evenly-spaced versus unevenly-spaced data, I can see no advantage to unevenly-spaced data as such, as far as accuracy of results is concerned. It is just that certain methods of data collection will naturally tend to produce uneven spacing in the data intervals, so I hold that it is prudent to provide the appropriate tools for handling such data where higher order polynomial approximation is called for. Naturally there will be a larger amount of computation required to handle the more complex expressions involved in uneven spacing. (I'm afraid I didn't go into this last point adequately in this thread; I just let the expressions speak for themselves.)

  I hope the above discussion was sufficiently vague to get me by in answering your question, Bruno.

Roger Stafford