Accuracy of integral function.

Hi

I am solving a long equation with q,k,phi,D are constants and x variable. I am using the the matlab keyword "integral" to integrate the equation. when I vary D from 10^-3 to 10^-5, I am supposed to get the final answer rho_ap vary from ~1 to ~0. However in my matlab code it gived my a fixed value of 1 when ever i vary D. Is the way i calculate my integration accurate. How accurate is integral function in matlab.

Thanks

      clear all
  clc
  format long

% definition of constants
rho0=3*(10^11);    
d= 1*10^-14;       
q=1.6*10^(-19);
k=25*10^-3;
Lb= 44*10^-6;
phi=90*10^-3;

f = @(x) (4*k/q)*( atanh(tanh(q*phi/(4*k))*exp(-x/Lb)) + atanh(tanh(q*phi/(4*k))*exp(-(d-x)/Lb)));

fun=@(x) 1./(cosh(q.*f(x)./k));

rho_ap= (2/d)*integral(fun,0,d/2);

Answer 1

[John D'Errico]

Sometimes we need to apply common sense. Is this the accuracy of the integral function, or something else?

WHENEVER you try to do an integration, optimization, anything, PLOT what you have. LOOK AT IT. Apply critical thought to what you have done so far. Just common sense. Don't just throw it into integral, and then wonder about the result.

You are integrating from 0 to d/2. So plot it over that domain.

ezplot(fun,[0 d/2])

Your function is essentially constant at 1 over that domain. Integral has no problem integrating the number 1 over any domain that you wish. This has absolutely nothing to do with integral, or the accuracy thereof.

In fact, if you look very carefully at the result of fun, it is IDENTICAL to 1 over that entire domain, down to the least significant bit of the result.

So you need to sit down and look more carefully at your function as you wrote it. If the computation is indeed correct, then you need to consider if what you are trying to do makes sense. If everything checks out, then you need to consider if it is possible to evaluate what you have written in double precision arithmetic. But if it is not possible to do so, then you need to rethink what you are doing.

Again. Common sense. NEVER just throw something at a utility without thinking about what you did, and if what you did makes any sense at all. Always plot anything that you can. If you can't plot it, then find some way to look at it. Does it make sense? Computers don't think for you, at least not yet. And once they start doing that, we may all be in trouble. Computers do exactly as you tell them to do, regardless of whether it is nonsensical.

Edit: here is a hint: LOOOK CAREFULLY at your constants. Think about them.

tan(q*phi/4*k)

What is the argument to the tan function?

q*phi/4*k
ans =
        9e-23

Yes, that small of a number. Of course, then you multiply it by

exp(-x/Lb)

where x is no larger than d/2.

exp(-(d/2)/Lb)
ans =
    1.161e-05

So the product is even smaller.

All of this is highly likely to produce meaningless garbage when computed in double precision arithmetic. Is it even correct? That I cannot know.

[friet]

Thanks for you comment. I am expected to see a little change of the final answer i.e., from 0.5-1. Si I apply "format long". Since the change in the final answer is very small, that why I ask about accuracy. /Thanks

[John D'Errico]

Again, this is not a question of the accuracy of integral. Fix the problem that Roger found. Then plot it. Verify that something is happening, as if it is not changing from 1, then you are wasting your time using double precision arithmetic.

Note that the difference is not that significant.

q*phi/(4*k)
ans =
     1.44e-19

Still a tiny number. For x==d/2,

tan(q*phi/(4*k))*exp(-(d/2)/Lb)
ans =
   1.6719e-24

What do you think atan of that number is?

Again, even if you have done nothing else wrong, you cannot solve this problem in double precision arithmetic. It still has nothing to do with integral.

[friet]

Hi John! thanks. I had a typo in the constant. Now i get it fix!! Thanks for your feedback and hint to check.

Answer 2

[Roger Stafford]

I see an error. The line

   f = @(x) (4*k/q)*( atan(tan(q*phi/4*k)*exp(-x/Lb))+ atan(tan(q*phi/4*k)*exp(-(d-x)/Lb)))

should be:

   f = @(x) (4*k/q)*( atan(tan(q*phi/(4*k))*exp(-x/Lb))+ atan(tan(q*phi/(4*k))*exp(-(d-x)/Lb)))

You were multiplying by k rather than dividing by it.

[friet]

Thanks Roger! i edit my question based on your suggestion. I also mis-typed tanh as tan and edit it also, however it didn't help:(

[friet]

Hi Roger! thanks. I had a typo in the constant. Now i get it fix!! Thanks for your feedback and hint to check.