Asked by Harsha K
on 13 Sep 2017

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!

Answer by John D'Errico
on 13 Sep 2017

Edited by John D'Errico
on 13 Sep 2017

Accepted Answer

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.

Sign in to comment.

Answer by 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.

Sign in to comment.

Answer by Harsha K
on 13 Sep 2017

Dear Teja, Dear John,

Thank you very much. I solved it in the following way:

I defined a function for finding it using symbolic variables and functions associated.

function h_opt = zero_area_optimum_d_finder(b_int, pd_def, R_def)

syms x y A d_opt;

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

A = int(y, x, -b_int/2, b_int/2);

h_opt = double(solve(A==0,d_opt));

end

Use the h_opt each time to shift my function/curve in the y-direction :-).

Please advice if there is an issue with the above implementation.

John D'Errico
on 13 Sep 2017

Harsha K
on 6 Jun 2019

Thank you very much John for the effort taken in writing that exhaustive and informative answer.

Just to requote your prediction and answer: In the other pair of regions, the integration problem will fail, and no real value of d_opt will exist.

This is exactly what happened. And I have caught this exception in my code and implemented a work around.

Appreciate a lot.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.