quadgk is wrong???

7 ビュー (過去 30 日間)
qingwei wang
qingwei wang 2019 年 3 月 18 日
編集済み: John D'Errico 2019 年 3 月 18 日
Wo=10, T=500, p=100, f(x) follows the normal distribution N (350,100);
1552897246(1).png
This is a function image, obviously an increase function, but it will be reduced to 0.1.png
This is the code
p=100;
wo=10;
T=500;
f1=@(x,mu)(x.*p.*normpdf(x,350,10));
Ewr=@(mu) quadgk(@(x) f1(x,mu),0,(mu-1/wo)*T)
B=0:1:1000;
Ewr_mu=arrayfun(@(mu) Ewr(mu),B);
plot(B,Ewr_mu)
grid on
title('Image of expected profit')
xlabel('mu')
ylabel('E(wr)')
Anyone can help me??Thanks!!!

採用された回答

John D'Errico
John D'Errico 2019 年 3 月 18 日
編集済み: John D'Errico 2019 年 3 月 18 日
Your normal distribution is one that is too narrowly defined around the mean of 350. In fact, it is effectively zero towithin tolerances outside of 350 +/- 60, as that would be 6 sigma in each direction. So what happens is for most limits of integration, quadgk sees something that seems to be zero everywhere it looks. Therefore the integral must be zero.
That is, over the domain [0,1000], that kernel can look to be an essential Dirac delta, centered at 350.
This is a common mistake made by people. When I see someone claim that a numerical integration routine is incorrect, I first would look to see if they have a narrow normal distribution in the kernel.
For example,see what happens when I do this:
x = 0:100:1000;
normpdf(x,350,10)
ans =
3.9404e-268 7.6539e-138 5.5307e-51 1.4867e-07 1.4867e-07 5.5307e-51 7.6539e-138 3.9404e-268 0 0 0
So, only if I look in a narrow region do I see anything that is even slightly different from 0. What is the integral of a function that is 0 over any domain? Zero. You always need to be careful when numerically integrating a function that looks like a delta function. And normal distributions tend to look exactly like that when viewed over a too wide interval.
This is the code you wrote:
p=100;
wo=10;
T=500;
f1=@(x,mu)(x.*p.*normpdf(x,350,10));
Ewr=@(mu) quadgk(@(x) f1(x,mu),0,(mu-1/wo)*T)
B=0:1:1000;
Ewr_mu=arrayfun(@(mu) Ewr(mu),B);
plot(B,Ewr_mu)
So, f1 is not a function of mu, which enters into the problem only as the upper limit of integration. Essentially, you are varying B from 0:1000, which then appears as mu.
But now look closely at what you do. I'll just use a limited set of values here.
mu = 0:100:1000
mu =
0 100 200 300 400 500 600 700 800 900 1000
(mu-1/wo)*T
ans =
-50 49950 99950 149950 199950 249950 299950 349950 399950 449950 499950
Now, think about those numbers as the upper limits of integration, for a function that looks VERY MUCH like a Direc delta in context.
This is a completely expected result, that quadgk will fail here.
  1 件のコメント
qingwei wang
qingwei wang 2019 年 3 月 18 日
I got it, thank you very much! ! !

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by