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);
采纳的回答
John D'Errico
2017-3-1
编辑：John D'Errico
2017-3-1
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.
3 个评论
friet
2017-3-1
John D'Errico
2017-3-1
编辑：John D'Errico
2017-3-1
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
2017-3-1
更多回答（1 个）
Roger Stafford
2017-3-1
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
