How to integrate using while loop?

9 visualizaciones (últimos 30 días)
Andrejus
Andrejus el 21 de Oct. de 2014
Comentada: Torsten el 3 de En. de 2017
Hi all,
I have a function: PM =@(w) (5/16*3^2.*(0.5578^4./w.^5).*exp(-5/4.*(w./0.5578).^(-4)));
its integral: A=integral(PM, 0, 5);
I want to divide the area under the function into equal area pieces. Each beeing of a size= A/n, where n is any number desired. To achieve this I was trying to use while loop, but did not succed. Additionally, I need to keep trace at which values of w the condtion is met.
Thanks.

Respuesta aceptada

Mike Hosea
Mike Hosea el 21 de Oct. de 2014
Editada: Mike Hosea el 21 de Oct. de 2014
I'm assuming this is homework. Let me give you a couple of hints.
1. You can make a function like this
q = @(b)integral(PM,0,b);
Then, as you know,
A = q(5)
But you can make a useful function like this
g = @(b)q(b) - A/5
This function g is zero when integral(PM,0,b) is 1/5 of integral(PM,0,5). Then you might want to look at
g = @(b)q(b) - (2/5)*A.
That function is zero when integral(PM,0,b) = (2/5)*A. And so forth...
2. FZERO is a function that might interest you.
  3 comentarios
Mike Hosea
Mike Hosea el 22 de Oct. de 2014
Editada: Mike Hosea el 22 de Oct. de 2014
It works fine in a loop. I just did it and calculated 100 points for exp(-x.*x) over the interval [0,10]. You can use a while loop, but then you have to increment the index variable yourself. It is easier to use a for loop on k = 1:n-1. Then the k-th w value is given by
g = @(w)q(w) - k*(A/n)
w(k) = fzero(g,[1e-10,5]);
I had to use 1e-10 instead of 0 for a left bound because your integrand blows up there. If you want to get fancy, you can use the previously found w value for the lower bound given to fzero instead of 1e-10.
The only reason I can see for using a while loop in particular is if you were going to implement the root-finding algorithm yourself, say using bisection.
Andrejus
Andrejus el 23 de Oct. de 2014
Editada: Andrejus el 23 de Oct. de 2014
Yey it works just as I wanted. Thank you very much!

Iniciar sesión para comentar.

Más respuestas (3)

Roger Stafford
Roger Stafford el 21 de Oct. de 2014
Actually, Andrejus, if you make a change of variable from w to x according to:
w = x^(-1/4),
you will see that your problem is suddenly converted into a much simpler problem that requires neither numerical integration nor the use of 'fzero'. Ordinary calculus will do the job.
  3 comentarios
Roger Stafford
Roger Stafford el 22 de Oct. de 2014
If the symbolic toolbox's 'int' is applied to the given function, the result is a simple analytic function. In attempting to divide this up into n equal parts the transformation w = x^(-1/4) would almost surely be found to be useful, so one way or another that transformation is begging to be performed, even if the approach is purely by way of Matlab.
Mike Hosea
Mike Hosea el 22 de Oct. de 2014
Let me say at the outset that you're probably right that this is a calculus homework exercise and not a MATLAB programming homework exercise. As a numerical analyst whose use of MATLAB goes back to the FORTRAN version, I do marvel that the Symbolic Toolbox is considered "purely MATLAB" around here. Whether it is or isn't is semantics. When I say "MATLAB", however, it means numerical. I've engineered plenty of test cases for numerical algorithms starting from a solution and then working "backwards" to arrive at a test problem. It's a way of engineering problems that I have an exact answer to for comparison purposes.

Iniciar sesión para comentar.


Roger Stafford
Roger Stafford el 22 de Oct. de 2014
In order to show that your problem need not involve time-consuming iterative methods here is the complete solution in four lines:
k1 = 5/16*3^2*(0.5578)^4;
k2 = 5/4*(0.5578)^4;
A = k1/k2/4*exp(-k2/5^4); % The total integral (area)
w = (1/5^4-log((1:n)'/n)/k2).^(-1/4); % Upper edges of equal area intervals
The justification of this code comes 1) from using 'int' (or ordinary calculus) to establish that the integral from 0 to t of k1/x^5*exp(-k2/x^4) is simply
k1/k2/4*exp(-k2/t^4)
and 2) that the solution to choosing intervals that divide this into n equal parts, that is, the w(m) such that
m/n = (k1/k2/4*exp(-k2/w(m)^4))/(k1/k2/4*exp(-k2/5^4))
is
w(m) = (1/5^4-log(m/n)/k2).^(-1/4)
by using either 'solve' or ordinary algebra.
  1 comentario
Andrejus
Andrejus el 23 de Oct. de 2014
Thank you for you attempt to help me out. I am sure that your approach is just as good as Mike's but this time I could get his idea faster. Thank you once again!

Iniciar sesión para comentar.


Roger Stafford
Roger Stafford el 24 de Oct. de 2014
Strictly speaking, this is an answer to a question you haven't asked, Andrejus, but I thought you might be interested in the fact that your request "to divide the area under the function into equal area pieces" is almost the equivalent of asking to fill the space below the curve with random points uniformly distributed area-wise (in a statistical sense.) Because the inverse of the integral of your function is easy to solve for, this can be done with the following code, which is almost identical to the solution I gave earlier to your original request. Ignoring the y values, this is also the solution to generating a random variable, x, whose probability density function (pdf) is the integrand you gave (after being normalized by being divided by the total area.)
k1 = 5/16*3^2*(0.5578)^4; % = 0.27227425026348
k2 = 5/4*(0.5578)^4; % = 0.12101077789488
n = 5000;
x = (1/5^4-log(rand(n,1))/k2).^(-1/4); % <-- The random (x,y) points
y = k1./x.^5.*exp(-k2./x.^4).*rand(n,1);
X = linspace(1e-10,5,500); % <-- The (X,Y) points in the curve above
Y = k1./X.^5.*exp(-k2./X.^4);
plot(X,Y,'r-',x,y,'y.')
  2 comentarios
Bambang Dewandaru
Bambang Dewandaru el 3 de En. de 2017
I tried to follow your step but even at this stage I did get my function work. Please help. function [ A, PM] = bd(k) a=1; PM= @(w) (a/pi*sin(w*pi/a)+w) A=@(k)integral(PM,0,k);
end
Torsten
Torsten el 3 de En. de 2017
A=integral(PM,0,k)
instead of
A=@(k)integral(PM,0,k)
Best wishes
Torsten.

Iniciar sesión para comentar.

Categorías

Más información sobre Logical en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by