Recursive method to get the differential not working.

Im trying to make a recursive method to get the n:th-order differential equation.

what i have currently is 2 methods im my .m file first one being the simple 1st order differential.

function func = differential(f) % callculates the n:th-order differential
    arguments
        f function_handle
    end
    h = 10^(-5);
    func = @(x)((f(x+h)-f(x))./h);
end

then im trying to use this in my recursive method

function output = differentialPower(f,n)
    arguments
        f function_handle
        n
    end
    if(n==0)
        output = f;
        return;
    else
        f = differentialPower(differential(f),n-1);
        output = f;
        return;
    end
end

to get the n:th differential

Problem i have is that my output will allways be either the original function (f)

or the first order differential of:

f = @(x)((f(x+h)-f(x))./h)

gooten from the differential method.

what i want to happen is that each time it gose deeper it will replace f(x) with ((f(x+h)-f(x))./h) and there for going deeper.

is this possible without using syms? or do i have to use the syms methods?

Answer 1

[Uday Pradhan]

Hi Hampus,

You will need to use the output function handle to evaluate the nth derivative approximation:

f = @(x) x^5;
df3 = differentialPower(f,3);
%Find approximation of f'''(1)
>> df3(1)

ans =

  60.174087934683492

This approximation will get highly inaccurate as you increase the order of the derivative. Instead, a neat way to calculate higher order derivatives is using the diff function.

[Hampus Toft]

Yes i noticed that the resault got worse and worse for each derivative when using a function like e^x

where dfn(0) always should be 0 but when using this approximation the answer started to giving answers ≠ 1...

I wanted to use this to calculate a Taylor Polynomial for e^(-x), but in the end i wrote a custom function to get the Taylor Polynomial only for e^(-x).

creating a function that would convert x^k/k! to a Pn(x) function for later use. Code for anyone intreasted:

% Method to give n:th Taylor Polynom String
function output = taylor(taylorFunction,n)
    fstring = "";
    for loop = 0:n
        if(loop==0)
            varfix = strrep(func2str(taylorFunction),'@(n,x)',"@(x)");
            fstring = strrep(varfix,'n',"(" + num2str(loop) + ")");
        else
            fstring = fstring + "+" + strrep(strrep(func2str(taylorFunction),'@(n,x)',""),'n',"(" + num2str(loop) + ")");
        end
    end
    output = str2func(fstring);
end

Usage would be:

% Taylor function for e^(-x)
fTaylor = @(n,x) ((-1).^n*(x.^n./factorial(n)));

% gets the Taylor Function of order n using the fTaylor
Pn = taylor(fTaylor,n)
PnFunction = stringFormat(Pn);

PnFunction is a optional poorly written function_handle -> latex text converter

function out = stringFormat(in)
    instring = func2str(in);
    
    displayString = strrep(instring,'@(x)',"");
    displayString = regexprep(displayString,"\.",""); % finds . and replace with nothing
    
    % removes redundant (-1)^n
    displayString = regexprep(displayString,"(\+*\(\(-1\)\^\((.?[02468])\)\*)","+");
    displayString = regexprep(displayString,"(\+*\(\(-1\)\^\((.?[13579])\)\*)","-");

    % formats exponetials to latex format i.e from 1^(15) to 1^{15} etc
    displayString = regexprep(displayString,"\^\((?<exponetial>.*?)\)","\^\{$<exponetial>}");

    % formats matlabs factorial(n) function into n!
    displayString = regexprep(displayString,"factorial\((?<factorial>.*?)\)","$<factorial>!");

    % an attempt to fix fractions in latex, proabably wrong but works in
    % this use case...
    expression = "\((?<numerator>.*?)\/\((?<denominator>.*?)\)([\+\-\*\!]?)\)\)";
    displayString = regexprep(displayString,expression,"\\frac{$<numerator>}{$<denominator>}");

    % dont remeber what this part even fixes :/
    out = regexprep(displayString,"^\+","");
end

The reason for trying to get n:th derivative numerical was due tou the assignment not allowing the usage os the syms toolbox, and i wanted a method to get the generic Taylor polynomial.

Anyway i'm going to consider this answard.