Problem [Matlab]

Prev: Relation between Hough & Radon Transform Matrices
Next: Help fixing Gosper curve (image not as it should be)

From: Roger Stafford on 4 May 2010 20:30

Beaver <marden.nomm(a)gmail.com> wrote in message <1606833116.71939.1273004986245.JavaMail.root(a)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

From: Beaver on 5 May 2010 14:01

Thanks to all. I find good advices.

Is this formula have some notes how this works?

From: Roger Stafford on 5 May 2010 21:53

Beaver <marden.nomm(a)gmail.com> wrote in message <1317322228.78895.1273096933766.JavaMail.root(a)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

From: Roger Stafford on 6 May 2010 04:15

"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <hrt7e0$drs$1(a)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

From: TideMan on 6 May 2010 04:24

On May 6, 8:15 pm, "Roger Stafford"
<ellieandrogerxy...(a)mindspring.com.invalid> wrote:
> "Roger Stafford" <ellieandrogerxy...(a)mindspring.com.invalid> wrote in message <hrt7e0$dr...(a)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.

First | Prev | Next | Last
Pages: 1 2 3 4
Prev: Relation between Hough & Radon Transform Matrices
Next: Help fixing Gosper curve (image not as it should be)