## How I do evaluate a function handle in other function handle

### capj (view profile)

on 15 Nov 2019 at 8:36
Latest activity Edited by Matt J

### Matt J (view profile)

on 9 Dec 2019 at 14:10
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

### Matt J (view profile)

on 15 Nov 2019 at 20:26
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.

R2019b

### Matt J (view profile)

on 15 Nov 2019 at 10:18

I think this line needs to be
g=@(x,y) f(x,y)*p_y(f,x,y) + p_x(f,x,y);

capj

### capj (view profile)

on 15 Nov 2019 at 17:52
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

### Matt J (view profile)

on 15 Nov 2019 at 17:58
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

### capj (view profile)

on 15 Nov 2019 at 19:01 where y=0.5 and x=0

### Steven Lord (view profile)

on 15 Nov 2019 at 14:37

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

### capj (view profile)

on 15 Nov 2019 at 15:04
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.

### Guillaume (view profile)

on 15 Nov 2019 at 16:21

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

### capj (view profile)

on 15 Nov 2019 at 16:32
Yes, I made a mistake with the y_p. How do you suggest the delta change would be, more specifically?
Guillaume

### Guillaume (view profile)

on 15 Nov 2019 at 16:49
x * 1e-4 or something similar would make more sense to me. That way the delta changes with the magnitude of x.
capj

### capj (view profile)

on 15 Nov 2019 at 17:32
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?

### Matt J (view profile)

on 15 Nov 2019 at 21:27

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.

Show 1 older comment
Matt J

### Matt J (view profile)

on 8 Dec 2019 at 2:55
Cubic splines have only 3 derivatives, higher order splines have more.
capj

### capj (view profile)

on 8 Dec 2019 at 3:05
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

### Matt J (view profile)

on 9 Dec 2019 at 13:44
I've never used it, but this File Exchange submission advertises spline fitting of any order
It also appears to contain functions for evaluating spline derivatives.