Moving a curve in one axis to ascertain area under the curve is zero

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 1

[John D'Errico]

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.

Answer 2

[Teja Muppirala]

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.

Answer 3

[Harsha K]

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]

You apparently did not bother to read the answer I spent a great deal of time writing. It explains exactly when you will have a problem.

[Harsha K]

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.