How I do evaluate a function handle in other function handle

I define the next function:

function [v] = f_deriv(f,x,y,n)
%f is a function of two variables
%x,y is a number
%n is the number of proccess that I need to do

p_y=@(f,x,y) (f(x,y+0.0001) - f(x,y-0.0001))/(2*0.0001);
p_x=@(f,x,y) (f(x+0.0001,y) - f(x-0.0001,y))/(2*0.0001);
y_p=f(x,y);


for i=1:n
    g=@(x,y) f(x,y)*p_y(f,x,y) + p_x(f,x,y);    
    f=@(x,y) g(x,y)
    
end

v=f(x,y)

end

This function has a input f, like (f=@((x,y) y -x^2 +1), and want to has as output $f^n(x,y)$ (the nth derivative of f, where y depends of x)

The steps that I follow to solve this problem are:

First, I define the function p_y and p_x that is the partial derivate aproximation of f. ()

Next, I define a y_p that is the function f evaluated in (x,y)

And, I define the function g, which is the (p_x + p_y*y_p)

Finally, I want to evaluate n times the function g_k in the same function, something like this:

or something like that:

g_1=@(x,y) f(x,y)*p_y(f,x,y) + p_x(f,x,y);    
g_2=@(x,y) g_1(x,y)*p_y(g_1,x,y) + p_x(g_1,x,y); 

.
.
.
g_n=@(x,y) g_n-1(x,y)*p_y(g_(n-1),x,y) + p_x(g_(n-1),x,y); 

But, when I use the code first I have incorrect results for n> 2. I know the algorithm is correct, because when I manually do the process, the results are correct

A example to try the code, is using f_deriv(@((x,y) y -x^2 +1,0,0.5,n), and for n=1, the result that would be correct is 1.5. For n=2,the result that would be correct is -0.5, for n=3. the result that would be correct is -0.5....

Someone could help me by seeing what the errors are in my code?, or does anyone know any better way to calculate partial derivatives without using the symbolic?

[Matt J]

I know the algorithm is correct, because when I manually do the process, the results are correct.

The algorithm is correct analytically, but is not stable numerically.

Answer 1

[Matt J]

I think this line needs to be

g=@(x,y) f(x,y)*p_y(f,x,y) + p_x(f,x,y); 

[Matt J]

Note that coding will be less bug-prone if you avoid reusing the symbols x and y for unrelated things.

[capj]

Thank you for your answer.

I think the suggestion is correct, but I get incorrect results for n = 3.(the result that would be correct is -0.5, and I get -4.5)

[Matt J]

I suggest you show the test you did of the code, so it can be reproduced.

[capj]

Sorry,

function [v] = f_dora(f,x0,y0,n)

p_y=@(f,x,y) (f(x,y+0.0001) - f(x,y-0.0001))/(2*0.0001);
p_x=@(f,x,y) (f(x+0.0001,y) - f(x-0.0001,y))/(2*0.0001);



for i=1:n
    g=@(x,y) f(x,y)*p_y(f,x,y) + p_x(f,x,y)
    f=@(x,y) g(x,y)
end
v=f(x0,y0);
end


 f_dora(@(x,y)y-x^2, 0 , 0.5 , 3)
 % -4.5, but the result that would be correct is -0.5

[Matt J]

I think I see a problem in the fundamental calculus. The derivative of f(x,y) with respect to x is

So, you have to know , which isn't made available in your code.

[capj]

In my case, I need that dy/dx = f(x,y). The code above is only for obtain all possible n-ths derived from f. Next, what I will proceed to do is that with these derivatives I calculate an approximate solution of the differential equation using the Taylor method.

[Matt J]

OK, but how do we see that the correct result at 0.,0.5, 3 is -0.5? Do we know the closed form of y(x) somehow?

[capj]

where y=0.5 and x=0

Answer 2

[Steven Lord]

g_1=@(x,y) y_p*p_y(f,x,y) + p_x(f,x,y);    
g_2=@(x,y) y_p*p_y(g_1,x,y) + p_x(g_1,x,y); 
.
.
.
g_n=@(x,y) y_p*p_y(g_(n-1),x,y) + p_x(g_(n-1),x,y); 
g_n(x,y)

Rather than using numbered variables, store your function handles in a cell array. That will make it easier for you to refer to previous functions.

f = @(x) x;
g = cell(1, 10);
g{1} = f;
for k = 2:10
    % Don't forget to _evaluate_ the previous function g{k-1} at x
    %
    % g{k} = @(x) k + g{k-1} WILL NOT work
    g{k} = @(x) k + g{k-1}(x);
end

g{10}([1 2 3]) % Effectively this returns [1 2 3] + sum(2:10)

[capj]

Using your suggestion I wrote the following code,

function [v] = f_doras(f,n)



g = cell(1,n);
g{1} = f
for k=2:n
    g{k}=@(x,y) g{k-1}(x,y)*(g{k-1}(x,y+0.001) - g{k-1}(x,y-0.001))/(2*0.001) + (g{k-1}(x+0.001,y) - g{k-1}(x-0.001,y))/(2*0.001);    
    
     
end
v=g{n};
end
a=f_dora(f,4)
a(0,0.5)
% -4.5, but the result that would be correct is -0.5

For n=2 and n=3, the results are correct, but for n=4 it's not right.

Answer 3

[Guillaume]

Your y_p never changes in your function. It's always the original . Shouldn't y_p be reevaluated at each step?

Your 0.0001 delta is iffy as well. It's ok when you're evaluating the derivative for but for small x (and y) it's not going to work so well. Should the delta be based on the magnitude of x (and y)?

[capj]

Thanks for answer,

Yes, I made a mistake with the y_p. How do you suggest the delta change would be, more specifically?

[Guillaume]

x * 1e-4 or something similar would make more sense to me. That way the delta changes with the magnitude of x.

[capj]

I don't understand your observation very well, for example, seeing the following code (and doing what I think I understood you)

function [v] = f_doras(f,n)

g = cell(1,n);
g{1} = f
for k=2:n
g{k}=@(x,y) (g{k-1}(x,y)*(g{k-1}(x,y+0.0001*y) - g{k-1}(x,y-0.0001*y)) + (g{k-1}(x+0.0001,y) - g{k-1}(x-0.0001,y)))/(2*0.0001*y);    
    
    
     
end
v=g{n};
end

I cannot make the magnitude change, because for example, if x = 0, then x * 0.0001 would be 0. If I do it with y, I think I have worse results than without making the change you suggested.

for example, using a= f_doras(@(x,y) y - x^2 + 1,3), and a(0,0.5) the result that I obtain is worse than using constant delta. Could you explain me more specifically your suggestion?

Answer 4

[Matt J]

What you probably have to do is use ode45 or similar to solve the differential equation dy/dx = f(x,y) on some interval. Then, fit f(x,y(x)) with an n-th order smoothing spline. Then, take the derivative of the spline at the point of interest.

[capj]

Thanks for the idea, but splines only have 3 derivaties.

[Matt J]

Cubic splines have only 3 derivatives, higher order splines have more.

[capj]

Is there an implementation in matlab of splines in a higher order,I have not found it, if it exists, you have given me the solution to my problem!

[Matt J]

I've never used it, but this File Exchange submission advertises spline fitting of any order

https://www.mathworks.com/matlabcentral/fileexchange/71225-splinefit

It also appears to contain functions for evaluating spline derivatives.