2 views (last 30 days)

Show older comments

Hi Folks,

Following is my function:

% where pd_def = defined constant which keeps changing for different situations

% b = pre-defined constant

% R = is calculated on the basis of pd_def

% x = x-axis variable which lets say assumes values in an interval [-b b]

% d_opt = variable which needs to take a value for defined parameters such that area under the

% curve 'y' is zero. This value once computed will be used to move the curve in y-direction to ensure that area under the

curve is zero.

R = b^2/(8*pd_def)+(pd_def/2);

y = R-sqrt(R^2-x^2) - (pd_def - d_opt);

I am having issues implementing this. Do I have to use optimisation toolbox to which I have no subscription. I would prefer an elegant solution without using the optimisation toolbox.

Please advice.

Cheers!

John D'Errico
on 13 Sep 2017

Edited: John D'Errico
on 13 Sep 2017

This answer gets complicated, because it looks forward to a problem that you have not foreseen, but you will probably trip over at some point.

First, it was not obvious to me that x is given already as a vector of values that span the interval [-b,b], or that x is essentially a symbolic variable, that just lives on the interval [-b,b].

In the first case, the integral will be done using a tool like trapz. In the latter case, integral works. Sort of. Read on...

Having said that, there is still a serious problem that I foresee, one that applies in either case.

If x lives on the interval [-b,b], then the integrand must exist on that interval. Depending on the values of b and pd_def, it is quite easy to have combinations of those variables such that sqrt(R^2-x^2) is a complex number on at least some pat of the interval [-b,b]. As I said, this is a serious problem.

The point is, that the integrand is only valid here when x lives on the interval [-R,R]. So, if R is less than b, the integration will fail. If R==b, then the integration has a singularity at each end of the interval, something mainly important here if a tool like trapz is employed.

So, what matters is when will it be true that R<b? The dividing line occurs at the point where R==b. Taken as a function of the unknown pd_def, we formulate the equation, that is quadratic in b, as a function of pd_def.

b^2/(8*pd_def)+(pd_def/2) == b

If we multiply by 8*pd_def, things become more clear.

b^2 + 4*pd_def^2 = 8*pd_def*b

As a function of pd_def, this will be a conic form, one that looks at first glance to be elliptical, but some further algebra is needed. Perhaps if I had used variables like x and y here, it would be more obvious. But x and y were already taken in the early parts of the problem statement.

So, next, back to that bit of algebra I said was needed. Is that the form of an ellipse, a hyperbola, what? Under what circumstances is there a problem? Some algebraic rearrangement is needed. We can re-write the equation for b and pd_def as

4*(pd_def - b)^2 - 3*b^2 == 0

You can verify that fact easily enough. TRY IT!

How does that help us? This is a difference of two perfect squares! So we can re-write it now as

(2*pd_def - 2*b - sqrt(3)*b))*(2*pd_def - 2*b + sqrt(3)*b)) == 0

How does that help us? It gives us a pair of intersecting straight lines. They intersect at the origin of course. but between the lines

2*pd_def - 2*b - sqrt(3)*b == 0

2*pd_def - 2*b + sqrt(3)*b == 0

Again, these are straight lines. Think of them as forming an x that crosses at the origin. They form 4 regions in the plane that extend out to infinity in all directions. In two of those regions, your integration will perform with no problem, so it has a solution. In the other pair of regions, the integration problem will fail, and no real value of d_opt will exist.

Of course, I cannot know what values your variables (b and pd_def) will take on. But Murphy's law tells me that you will have a problem along the line. And then you will come back here, frantic, asking why your simple solution fails. What did you do wrong? In that case, carefully re-read my answer.

Teja Muppirala
on 13 Sep 2017

Calculate the value of the integral when d_opt = 0, and then use that to set d_opt to cancel the integral out.

b = 0.1;

pd_def = 1;

R = b^2/(8*pd_def)+(pd_def/2);

y = @(x) R-sqrt(R^2-x.^2) - (pd_def); % First set d_opt = zero

val = integral(y,-b,b)

"val" contains the value of the integral when d_opt is zero.

Then set d_opt = -val/(2*b) , where 2*b is the length of the integral.

d_opt = -val/(2*b)

y = @(x) R-sqrt(R^2-x.^2) - (pd_def - d_opt);

areaUnderY = integral(y,-b,b)

This gives d_opt =

0.996654840715272

areaUnderY =

-4.838252194716564e-18

Basically zero.

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

Start Hunting!