Home › Forums › Get Help › How do I ...? with Maple

Impossible to solve an integral?

Posted on 2008-01-08 07:29 By Xan Guillen (0)
Categories:

Hi, i have been fighting with Maple 11 for a long time and i would like to know if i am doing something wrong or the integral is so complex for the program... I really need to solve this problem soon if i would like to have a successfully Master Thesis :)

During the degree i study Maple but i forgot a lot :( Here is the issue:

I am working in telecommunication issues and calculating the output power i get this integral:

P[through] := int(2.000000000*10^8*(.51-1.428285686*10^(-0.1600000000e-1*Pi)*sqrt(1-k[2]^2)*cos(0.2473600e-2*Pi^2/lambda)+(10^(-0.1600000000e-1*Pi))^2*(1-k[2]^2))/((1-.7141428429*10^(-0.1600000000e-1*Pi)*sqrt(1-k[2]^2))^2*(1+2.856571372*10^(-0.1600000000e-1*Pi)*sqrt(1-k[2]^2)*sin(0.1236800e-2*Pi^2/lambda)^2/(1-.7141428429*10^(-0.1600000000e-1*Pi)*sqrt(1-k[2]^2))^2)), lambda = 6.30*10^(-7) .. 6.35*10^(-7))

There is only still one variable called k[2], and i would like to know the value of k[2] that maximize the P[through], being k[2] a real value between 0 and 1 (i think i have tried everything with assuming...)

Is there a way to do this and also to see the expression of P(k[2]) ? Thanks for your help!

I have also another problem more strange...

When i plot a graph in a range of 50nm the graph is not correct for certain ranges, so i plot again the same function in that range only (smaller range this time), but it is correct in this second graph! I have tried to put higher Digits, but it does not affect anything... Some help? Could be this affecting the integral in this interval?

‹ Looking for: Laplace lectures using Maple Simple for loop ›
John Fredsted's picture

Integral

Comment on 2008-01-08 14:06 from John Fredsted (Maple Rating 4 566)

Unless I am mistaken, your integral is proportional to the definite version of

int((a + b*cos(c/lambda)) / (1 + d*sin(e/lambda)^2),lambda);

where a, b, c, d, and e are constants (depending on k[2]). On face value, Maple is unable to integrate it, giving just

<br /> int((a+b*cos(c/lambda))/(1+d*sin(e/lambda)^2), lambda)<br />

But maybe some hardcore mathematicians (I am only a physicist) can help you out here.

Robert Israel's picture

integral

Comment on 2008-01-08 14:07 from Robert Israel (Maple Rating 4 1555)

The integral
int((a+b*cos(c/lambda))/(1+d*sin(e/lambda)^2), lambda) has almost certainly no closed-form solution in general. It looks to me in Xan's case c = 2*e. I think there may be relations between the other constants as well, e.g.
1.428285686 is presumably meant to be 2 * .7141428429. If you want to get closed-form solutions, there's a better chance if you leave the constants in symbolic form rather than using floating point. Not that it's very likely that there's a closed form solution. Still, if you could do reliable numerical evaluations of the integral you should be able to use numerical methods to find the maximum. But the integrand is oscillating pretty wildly in this interval, so numerical evaluation is going to be hard. A higher setting of Digits is almost certainly needed; in fact, it may already be too late as the result could be quite sensitive to the roundoff error that has already occurred.

Robert Israel's picture

Additional thoughts

Comment on 2008-01-08 15:35 from Robert Israel (Maple Rating 4 1555)

After using sin(e/lambda)^2 = 1/2 - 1/2*cos(c/lambda) and substituting sqrt(1-k[2]^2) = r, I find that the integral in question is a constant multiple of
Int((1-alpha*r*cos(c/lambda)+beta*r^2)/(1+gamma*r^2-delta*r*cos(c/lambda)),lambda )
with alpha = 2.494476546, beta = 1.555603312, gamma = .4046124214,
delta = 1.272183039. The change of variables
c/lambda=t makes the definite integral into
a constant multiple of
Int((1-alpha*r*cos(t)+beta*r^2)/(1+gamma*r^2-delta*r*cos(t))/t^2,t = 38446.38339 .. 38751.51341)
Note that the endpoints are large and relatively close together (though their difference is nearly 100 pi). So over this interval the factor 1/t^2 is nearly constant, and maximizing the integral should be mainly a matter of maximizing
S = Int((1-alpha*r*cos(t)+beta*r^2)/(1+gamma*r^2-delta*r*cos(t)),t=0..2*Pi)
Now
Int(1/(1-k*cos(t)),t=0..2*Pi) = 2*Pi/sqrt(1-k^2) for abs(k) < 1
so I get S = -2*Pi*(-alpha*(1+2*gamma*r^2-delta^2*r^2+gamma^2*r^4)^(1/2)+alpha+alpha*gamma*r^2-delta-delta*beta*r^2)/delta/(1+2*gamma*r^2-delta^2*r^2+gamma^2*r^4)^(1/2)<br />
Differentiating that with respect to r, I get r=0 and (+/-)1/(alpha*gamma*delta-delta^2*beta+2*beta*gamma-2*gamma^2)*((alpha*gamma*delta-delta^2*beta+2*beta*gamma-2*gamma^2)*(-delta^2-2*beta+alpha*delta+2*gamma))^(1/2)
which would be outside the interval [0,1]. So I suspect that the maximum of P[through] is either at k[2] = 0 or 1.

Numeric solution?

Comment on 2008-01-08 16:50 from DJKeenan (Maple Rating 2 191)

I agree with Robert Israel, that it might be too late to do numeric integration, because round-off error might already have occurred (in obtaining the constants). Having said that, can we try the following?

f:= unapply(P[through], k[2]);
plot(f, 0 .. 1);

which indicates that the function is decreasing on the interval, and thus the maximum is indeed at k[2]=0. You can test that via

Optimization:-Maximize(f, 0 .. 1, method='branchandbound');

Regarding increasing Digits, this does not affect the plot command, which tends to use hardware floats for evaluation. (There are ways around this, via defining a function with its own Digits.)

Re: Numeric solution?

Comment on 2008-01-09 03:59 from pantole (Maple Rating 1 34)

If it is a matter of evaluating numerically the integral, then you may test it using GlobalOptimization package (using branchandbound option), as DJKeenan suggested:

> restart:
> integral:=proc(k)
if type (k, numeric) then
int(2.000000000*10^8*(.51-1.428285686*10^(-0.1600000000e-1*Pi)*sqrt(1-k^2)*cos(0.2473600e-2*Pi^2/lambda)+(10^(-0.1600000000e-1*Pi))^2*(1-k^2))/((1-.7141428429*10^(-0.1600000000e-1*Pi)*sqrt(1-k^2))^2*(1+2.856571372*10^(-0.1600000000e-1*Pi)*sqrt(1-k^2)*sin(0.1236800e-2*Pi^2/lambda)^2/(1-.7141428429*10^(-0.1600000000e-1*Pi)*sqrt(1-k^2))^2)), lambda = 6.30*10^(-7) .. 6.35*10^(-7))
else 'integral'(k) end if; end proc:
> with(GlobalOptimization): infolevel[GlobalOptimization]:=3:
> GlobalSolve(integral(k),k=0..1, maximize = true, timelimit=57600, method=branchandbound);

The answer returned after 25 minutes (in a 1 GHz processor) is:

GlobalSolve: calling NLP solver
SolveGeneral: calling global optimization solver
SolveGeneral: number of problem variables 1
SolveGeneral: number of nonlinear inequality constraints 0
SolveGeneral: number of nonlinear equality constraints 0
SolveGeneral: method branchandbound
SolveGeneral: merit function evaluation limit 1000
SolveGeneral: non-improving merit function evaluation limit 200
SolveGeneral: constraint penalty multiplier 100.0
SolveGeneral: target merit function value 0.10e11
SolveGeneral: local search target objective function value 0.10e11
SolveGeneral: local search feasibility tolerance 0.10e-5
SolveGeneral: local search optimality tolerance 0.10e-5
SolveGeneral: time limit in seconds 57600
SolveGeneral: trying evalhf mode
SolveGeneral: trying evalf mode
LGORUN: total number of function evaluations 295
LGORUN: runtime in external solver 1583.
LGORUN: maximum constraint infeasibility 0.

[0.829292352751139994, [k = 0.]]

That is, GlobalSolve returns the lower bound. If theoretical considerations must be included, then you should refer to Robert's answer and compare P(0)/P(1) [having P function of k] or something like that.

I understand that few people have the Global Optimization toolbox, and it is rather inconvenient to post an example with that. However, I am posting the code because I believe that using the add-on package [Optimization] is not so straightforward as the above.

Hope that Maplesoft will consider the possibility of integrating GOT in the forthcoming versions of Maple.

pantole

acer's picture

Optimization

Comment on 2008-01-09 04:32 from acer (Maple Rating 4 1410)

Actually, DJKeenan suggested something like,

Optimization:-Maximize(integral,0..1,method='branchandbound');

which also returns [0.829292352740012, [0.]].

The purely local solver invoked as Optimization:-Maximize(integral,0..1) also gave that result for me, even though the function may not be monotonic (if the apparent behaviour between 0.99999 and 1.0 is accurate).

acer

Agreed

Comment on 2008-01-09 04:43 from pantole (Maple Rating 1 34)

Thanks for the correction. The built-in [Optimization] package is indeed more straightforward than the GOT, which I didn't test. What I had in mind is a previous post about an integral, which I briefly tested using the Optimization package.

Now I have to reconsider also my previous entry...

pantole

Axel Vogt's picture

r ~ 0.0024

Comment on 2008-01-09 17:45 from Axel Vogt (Maple Rating 4 994)

Following Robert Israel's way I use

L:=[alpha = 2.494476546, beta = 1.555603312, gamma = .4046124214, delta = 1.272183039];
K := (1-alpha*r*cos(t)+beta*r^2)/(1+gamma*r^2-delta*r*cos(t))/t^2; # his new integrand

Making the bounds rational (still working with that ugly 10 Digits ...) now use

 Int(K,t = 610605461/15882 .. 3875151341/100000, method=_d01akc);
 eval(%,L);
 plot(%,r=5*1e-4 .. 5*1e-3, numpoints=400);

to guess from the graphics (and its symmetry) that the minimum is the middle of
0.003 and 0.0018

The different result may be caused by

'eval(K, t = t + 2*Pi)/K': '%'=%;

                     K|                    2
                      |t = t + 2 Pi       t
                     -------------- = -----------
                           K                    2
                                      (t + 2 Pi)

Now start over, slicing (and using better exactness) should fill some gap in the thesis.

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}