@mehdi jafari

Please excuse me @Preben Alsholm if I take the liberty of intervening in this discussion.

To @mehdi jafari :

The cpu time seems (to me) reasonable given the complexity of A1 and A2

CodeTools:-Usage( plot(w-> int(eval(A1+A2,omega=w),theta=0..2*Pi, numeric, digits=5), 0..50, color=red) );
memory used=314.85MiB, alloc change=0 bytes, cpu time=4.79s, real time=4.64s, gc time=224.68ms

Given the symmetries of A1+A2 one can reduce it this way

CodeTools:-Usage(plot(w-> int(eval(A1+A2,omega=w),theta=0..2*Pi, numeric, digits=5), 0..50)):
CodeTools:-Usage(plot(w->4*int(eval(A1+A2,omega=w),theta=0..Pi/2., numeric, digits=5), 0..50)):
memory used=328.38MiB, alloc change=36.00MiB, cpu time=5.13s, real time=5.08s, gc time=179.32ms
memory used=215.62MiB, alloc change=40.00MiB, cpu time=2.90s, real time=2.84s, gc time=131.88ms

As A1+A2 contains a "*boundary layer of length of order 1e-4*" located around omega=16 (the integral is numerically null for omega < 16 - 1e-4)...

taylor(op(2, A1)+op(2, A2), omega=16, 3);
plot(w->int(eval(A1+A2,omega=w),theta=0..2*Pi,numeric,digits=5),16-2e-5..16+2e-5);

... one may think that computing the integral for omega >= 16 (maybe this value can be adjusted) is enough. What I propose is to plot this.

plot(w->4*int(eval(A1+A2,omega=w),theta=0..Pi/2., numeric, digits=5), 16-4e-5..50)

To realize an objective (I think) comparison, one must evaluate each plot with the same number of points.

Doing

pp := plot(w-> int(eval(A1+A2,omega=w),theta=0..2*Pi, numeric, digits=5), 0..50):
plottools:-getdata(pp)[3]

shows 231 values of omega are used in the range [0..50].

Let's fix the omega step to the same value 50/230 and evaluate the integration time alone (plotting time discarded). One obtains:

J := proc() seq(int(eval(A1+A2,omega=w),theta=0..2*Pi, numeric, digits=5), w=0..50, 50/230): end proc:
CodeTools:-Usage(J(), output='cputime', iterations=10):
memory used=180.73MiB, alloc change=0 bytes, cpu time=2.24s, real time=2.13s, gc time=168.40ms
K := proc() seq(4*int(eval(A1+A2,omega=w),theta=0..Pi/2, numeric, digits=5), w=0..50, 50/230): end proc:
CodeTools:-Usage(K(), output='cputime', iterations=10):
memory used=162.65MiB, alloc change=0 bytes, cpu time=1.96s, real time=1.85s, gc time=152.76ms
Ka := proc() seq(4*int(eval(A1+A2,omega=w),theta=0..Pi/2, numeric, digits=5), w=16..50, 50/231): end proc:
CodeTools:-Usage(Ka(), output='cputime', iterations=10):
memory used=57.89MiB, alloc change=0 bytes, cpu time=718.70ms, real time=675.30ms, gc time=61.45ms

Comparing the cpu times between J and Ka shows the latter is three times faster than the former.

Next refinement (at least for me): precalculate the list of points [x, int(...)] and plot this list.

Kb := proc()
local u:=[ seq([w, 4*int(eval(A1+A2,omega=w),theta=0..Pi/2, numeric, digits=5)], w=16..50, 50/231) ]:
plot(u):
end proc:
CodeTools:-Usage(Kb());
memory used=57.99MiB, alloc change=0 bytes, cpu time=879.00ms, real time=764.00ms, gc time=174.41ms

We are then 5 times faster than the original version given at the top ogfthis reply.

Probably better performances can be obtained by:

- choosing carefuly the integration method and its options
- adjusting the "omega step" (beyond, let's say, omega=30, the definite integral slighly evolves)

To end this: believe it would be more important to verify that the the results are not dependent of the integration method, of its options, and of the omega step.

For instance A1+A2 has discontinuities for all omega (at least for values right to the "boundary layer")r

discont(eval(A1+A2, omega=16.0001), theta);
{-0.5004339588e-1+3.141592654*_Z29, 0.5004339588e-1+3.141592654*_Z27, Pi*_Z67+(1/2)*Pi}

Are these discontinuities correctly handled with the integration procedure? Does this matter?

printf("%1.3e", evalf(eval(eval(A1+A2, omega=16.0001), theta=0.05004339586)))
1.600e+05

Interestingly the integral cannot be computed for **omega=20.099742 .. 20.099751** ; the peak obtained in the original plot is hugely underestimated:

int(eval(A1+A2,omega=20.099752),theta=0..Pi, numeric, digits=5)
2.3543