How to calculate derivative and then apply limit in matlab

How to calculate diff((x/(exp(x)-1)),x,n) and then apply the limit at x=0 Where n=11,12,13,14,15 I am getting upto 1 to 10 but not from 11 onwards... please anybody help me...

Answer 1

[Deepak Gupta]

Hi Ravikiran,

I have used a for loop.

syms f(x) x; 
f(x) = x/(exp(x)-1);
g = f;
limitg = sym(zeros(15, 1));
for n = 1:15
    g = diff(g);
    limitg(n) = limit(g, x, 0);
end

Code may throw error of "Devide by Zero".

Thanks,

Deepak

[Ravikiran Mundewadi]

I am getting error in command window.

[Ravikiran Mundewadi]

The given function is f(x)=x/(exp(x)-1); For first derivative we can do and so on upto 10th derivative i can do and apply the limit for every derivative at x=0. Next i am not getting the value for 11th derivative and applying the limit at x=0.

[Ameer Hamza]

Ravikiran, what do you mean by "not getting the value"? Are you getting an error?

[Ravikiran Mundewadi]

Ya while typing the above code i am getting error.

[Ravikiran Mundewadi]

While derivative of the given above function and apply the limit at x=0 will get the value -1/2. Similarly for second derivative and then apply the limit will get the value 1/6. Ans so on upto 10th derivative we can get, but for 11th, 12th, 13th,14th and so on will not get the value. So what to do?

[Ameer Hamza]

"will not get the value"??? I guess you are getting zero, right???

[Ravikiran Mundewadi]

Ya u r right... Y?

[Ravikiran Mundewadi]

Actually i have get the values for 1th derivative is zero, 12th derivative is -691÷2730, 13th derivative is zero, 14th derivative is 7/6, 15th derivative is zero, and so on...

[Ravikiran Mundewadi]

Actually we need to get these values...

[Ravikiran Mundewadi]

How it is possible in matlab?

[Ameer Hamza]

Ravikiran, How did you calculate these values? Did you calculate it in some other software?

[Ravikiran Mundewadi]

In Matlab

[Ameer Hamza]

I mean how do you know that for 12th derivative the value is -691/2730?

[Deepak Gupta]

Yes Ravikiran, How did you get the values for n = 11 to 15?

What i noticed is that for values greater than n=10, matlab fails to find a convergence for the limit x->0 of the expression and hence returns the expression itself. Which obviously can't be converted to double and therefore we observe the error.

Can you explain how did you get the values for n=11 to 15?

[Ameer Hamza]

Deepak, for n>10 MATLAB does not automatically give a simplified answer. But if you use simplify() instead of double() you will get a numeric answer.

syms x

y = x/(exp(x)-1);
D12y = simplify(limit(diff(y, 12), x, 0));

Result:

D12y =
0

Therefore, I am a bit confused. If Ravikiran claims that the correct answer is -691/2730, then this might be a bug in the symbolic engine.

[Deepak Gupta]

Ameer, Ravikiran's claim may be correct but not sure how's he getting to the result. Using simplify results will be different because simplify reduces the expression by removing less significant terms.

Using simplify, we can calculate results upto n=15 and ans is :

[ -1/2, 1/6, 0, -1/30, 0, 1/42, 0, 0, 0, 0, 0, 0, 0, 0, 0]

But without using simplify upto n = 10, ans is:

[ -1/2, 1/6, 0, -1/30, 0, 1/42, 0, -1/30, 0, 5/66]

So you can see results vary with and without using simplify.

In my understanding, For longer expressions, i.e. in this case expression for n=11 is

11242*exp(2*x))/(exp(x) - 1)^3 - (615780*exp(3*x))/(exp(x) - 1)^4 + (9003720*exp(4*x))/(exp(x) - 1)^5 - (56133000*exp(5*x))/(exp(x) - 1)^6 + (180789840*exp(6*x))/(exp(x) - 1)^7 - (325987200*exp(7*x))/(exp(x) - 1)^8 + (332640000*exp(8*x))/(exp(x) - 1)^9 - (179625600*exp(9*x))/(exp(x) - 1)^10 + (39916800*exp(10*x))/(exp(x) - 1)^11 - (11*exp(x))/(exp(x) - 1)^2 - (x*exp(x))/(exp(x) - 1)^2 + (2046*x*exp(2*x))/(exp(x) - 1)^3 - (171006*x*exp(3*x))/(exp(x) - 1)^4 + (3498000*x*exp(4*x))/(exp(x) - 1)^5 - (29607600*x*exp(5*x))/(exp(x) - 1)^6 + (129230640*x*exp(6*x))/(exp(x) - 1)^7 - (322494480*x*exp(7*x))/(exp(x) - 1)^8 + (479001600*x*exp(8*x))/(exp(x) - 1)^9 - (419126400*x*exp(9*x))/(exp(x) - 1)^10 + (199584000*x*exp(10*x))/(exp(x) - 1)^11 - (39916800*x*exp(11*x))/(exp(x) - 1)^12

Matlab fails to find convergence point and hence returns the expression itself.

[Ameer Hamza]

I get the same answer, with or without simplify(). I am using R2020a, and from recent questions on this forum, I have noticed that they have improved Symbolic toolbox in this release.

syms x
y = x/(exp(x)-1);

for i=1:15
    Dy = simplify(limit(diff(y, i), x, 0));
    if Dy == 0
        fprintf('i = %d,\t Dy = 0\n', i);
    else
        [n, d] = rat(Dy);
        fprintf('i = %d,\t Dy = %d/%d\n', i, n, d);
    end
end

Result:

i = 1,	 Dy = -1/2
i = 2,	 Dy = 1/6
i = 3,	 Dy = 0
i = 4,	 Dy = -1/30
i = 5,	 Dy = 0
i = 6,	 Dy = 1/42
i = 7,	 Dy = 0
i = 8,	 Dy = -1/30
i = 9,	 Dy = 0
i = 10,	 Dy = 5/66
i = 11,	 Dy = 0
i = 12,	 Dy = 0
i = 13,	 Dy = 0
i = 14,	 Dy = 0
i = 15,	 Dy = 0

[Ameer Hamza]

It turns out OP's claim is correct: https://www.wolframalpha.com/input/?i=limit+x-%3E0+d%5E12%2Fdx%5E12+x%2F%28e%5Ex-1%29

This is one of those situations where limitations of MATLAB symbolic toolbox become visible.

Answer 2

[Ameer Hamza]

As discussed in the comment to Deepak's answer. This is a limitation of of MATLAB's symbolic engine. MATLAB calculates the limit for n-th order derivatives for n>10 to be zero

syms x
y = x/(exp(x)-1);

for i=1:15
    Dy = simplify(limit(diff(y, i), x, 0));
    if Dy == 0
        fprintf('i = %d,\t Dy = 0\n', i);
    else
        [n, d] = rat(Dy);
        fprintf('i = %d,\t Dy = %d/%d\n', i, n, d);
    end
end

Result

i = 1,	 Dy = -1/2
i = 2,	 Dy = 1/6
i = 3,	 Dy = 0
i = 4,	 Dy = -1/30
i = 5,	 Dy = 0
i = 6,	 Dy = 1/42
i = 7,	 Dy = 0
i = 8,	 Dy = -1/30
i = 9,	 Dy = 0
i = 10,	 Dy = 5/66
i = 11,	 Dy = 0
i = 12,	 Dy = 0
i = 13,	 Dy = 0
i = 14,	 Dy = 0
i = 15,	 Dy = 0

But Wolfram Alpha is able to calculate the limit: https://www.wolframalpha.com/input/?i=limit+x-%3E0+d%5E12%2Fdx%5E12+x%2F%28e%5Ex-1%29. Even this FEX submission by John: https://www.mathworks.com/matlabcentral/fileexchange/20058-adaptive-numerical-limit-and-residue-estimation is not able to converge to a limit. I guess there is nothing much you can do about in MATLAB, other than writing your own closed-form solution if it is possible.

[Ravikiran Mundewadi]

Thanks to Deepak and Ameera. In wolffarm i tried to get the values, but upto 21sth derivative and 22nd and soon... its not displaying , its takes lot of time.

Actually its a concept of finding the Bernoulli numbers. whatever we tried its about the definition of Bernoulli number. i.e., Limit x tends to 0 d^n/dx^n (x/e^x-1).

[Ameer Hamza]

Ok. I got it now. This is the definition of Bernoulli number, but to generate it, you don't need to use this definition. There are several other efficient algorithms. Matlab has a built-in function to generate them. Try this

bernoulli(sym(0:30))

This require Symbolic toolbox.

[Ravikiran Mundewadi]

Oh superb... Thank you so much Ameer... I got it...

[Ameer Hamza]

I am glad to be of help.

[Ravikiran Mundewadi]

Be in touch... Thank you...