Matlab derivative approximation - can anyone help?

I've been puzzling over this last problem in my assignments for the past week or so, and still with no end in sight. Some of the folks here on MathWorks have given me a few pointers, but I still can't make this code work.

   function r=myderivative(f,a,tol)
   y = 0;
   maxIters = 500;
   n = 1;
   while n<=maxIters
       h = 1/n;
       while (1/n+1)-h<tol;
           y = (f(a+h)-f(a))/h;
       end
       n = n + 1;
       y = (f(a+h)-f(a))/h;
   end
   r=y;
   end

The problem in question instructed me to make an m-file mimicking a derivative function, which would apply (f(a+h)−f(a))/h to f(a) again and again (with h turning into 1, then 1/2, then 1/3...) until the successive values dipped under the "tol" value.

However, this current code skips over pretty much all the stuff after "y=0" and spits out 0 as the value for every set of (f,a,tol).

Could anyone spare a bit of help? I would greatly appreciate it.

[dpb]

Well, work thru the results of what you've written...for example, look at the following that displays the condition you've calculated for the comparison to tol in the second loop.

>> maxIters=3;
>> n=1;
>> while n<=maxIters
     h = 1/n;
     (1/n+1)-h
     n = n + 1;
   end

I just took out the computations and looked at the looping construct. You don't give a value for tol but note that your comparison value never changes as you've got it defined.

Answer 1

[Star Strider]

You need to track the successive values for the derivative y and compare them with the current estimate. The derivative calculation has to go before the second while block, and the second while block has to compare the previous value of the derivative with the current value. As long as that difference is >=tol, replace the previous value of the derivative with the current value and continue.

[Chong]

Okay, I'll try that. I'll Accept your answer as soon as I can get this to work. Wish me luck!

[Star Strider]

I actually got it to work with the changes I described, so I was confident in suggesting them. I’d have posted my solution, but I don’t want to deprive you of the satisfaction of getting it to work yourself! I’ll post it later if you like.

(This is relatively straightforward, and isn’t like the L’Hospital example, where you could not be expected to know that the anonymous functions lost their definitions as functions in the symbolic diff statement, and that you needed to use subs to make it work.)

[Chong]

Alright, I admit defeat.

I've been trying to take your suggestions into account for the past hour (shifting the derivative in front of the second while block, among other things), but the code still won't work properly. I'd like to post it here myself, but my Matlab seems to have frozen up entirely, so I can't even extract it.

Could you please simply post your solution?

[Star Strider]

OK. Here’s the ‘guts’ of the function:

     yp = 0;
     maxIters = 500;
     n = 1;
     while n<=maxIters
         h = 1/n;
         y = (f(a+h)-f(a))/h;
         while yp - y > tol;
             yp = y;
         end
         n = n + 1;
     end
     r=y;

You only needed to make a few changes.

[Chong]

Thank you for the effort, but it still looks like I'm missing something. I'm still not getting the "right" answers when I plug my sample data in.

For instance, the string (@(x) sin(x),pi/4,0.1) is supposed to produce a value of 0.577320817064314, but the code still gives me a value of 0.7064. I think the problem is that the code won't recognize my tolerance value for some reason (I get the exact same answer when I change the Tol to 0.01, or 0.001).

[Star Strider]

If I remember correctly, the derivative of sin(x) is cos(x), and cos(pi/4)=0.7071...

The code seems to be correct.

[Chong]

Ah, but the key here is that the code is meant to represent a rough approximation of the derivative, not a perfect one. Rough to the level of "let's do the y = (f(a+h)-f(a))/h by hand a couple dozen/hundred times and call it a day".

Thank you for all the help you've given me, but I think I'll try tackling this one on my own for a couple hours. I'll holler when I need more help.

Peace out.

[Star Strider]

It can be as rough as you want it to be, but has an appropriately low error at a tolerance of 1E-10 or so. For me, the whole idea of an iterative approximation scheme is to generate as accurate and precise an estimate as possible.

[Chong]

For the record, here are the directions verbatim:

Let's both work on it, and see who can get the right coding first. The acid test is still:

d = myderivative(@(x) sin(x),pi/4,0.1) = 0.577320817064314

May the best man win!

[Star Strider]

How did you come up with 0.577320817064314 as cos(pi/4)? That’s just never going to work because cos(pi/4)=0.7071...

If that’s the goal, I yield.

[Chong]

It's really not as complex as it looks.

The goal is basically to write a code that will apply [f(a+h)-f(a)]/h to any combination of (f(a),a,tol). The code is in a while loop, and will only stop when the difference between two successive h-values is smaller than the "tol" value.

For instance, I got 0.577320817064314 when I did [sin(pi/4+1/3) - sin(pi/4)]/(1/3). In this case, h is 1/3 because 1/2-1/3 is the last number that's bigger than the "tol" value, 0.1 (check for yourself: 1/3-1/4 is the first h-value that's smaller than 0.1).

Hope this helped.

(And no, I haven't gotten the code yet, either. Onward!)

[Image Analyst]

I don't understand your challenge. It's your homework. If you can do it, then do it. Why should Star waste his time doing something that you're already going to do by yourself anyway? I'm sure he could have done it already instead of giving you hints, if had wanted to.

[Chong]

Actually, we already turned in the thing a couple hours back. I didn't get this particular problem, but it was only a small piece of the homework.

Now I'm just trying to crack the thing for fun.

[Star Strider]

It actually appeared a week ago here, posted by a classmate. I solved it then for fun, essentially in the time it took to type it, but didn’t post it as an Answer to that Question. If the derivative of sin(pi/4) was supposed to be 0.57732..., it was a trick question, and no one could have gotten it.

[Chong]

Alright, after many hours, much brainstorming, and the sacrifice of a hobo, I finally got the code to work it the way I want. Preserved here for posterity:

function r = myderivative(f,a,tol)
h = 1;
oldapprox = (f(a+h)-f(a))/h;
h = h*(1/h)*(1/(1/h+1));
approx = (f(a+h)-f(a))/h;
while (abs(approx-oldapprox)>tol)
    oldapprox = approx;
    h=h*(1/h)*(1/(1/h+1));
    approx = (f(a+h)-f(a))/h;
end
r = approx;
end

(The real trick is figuring out how to get from 1/n to 1/(n+1) without introducing another bloody variable.)

[Star Strider]

My original code (from a week or so ago):

tol = 1E-12;
dfdx = @(f,a,h) (f(a+h)-f(a))./h;
dd = Inf;
dprv = 0;
k1 = 1;
while dd > tol
    h = 1/k1;
    df = dfdx(f,a,h);
    dd = df-dprv;
    dprv = df;
    k1 = k1+1;
end
df      % Result

produces:

df =

     707.1039e-003