## 5906 Reputation

19 years, 305 days
Munich, Bavaria, Germany

## for example...

Re(...)=0, Im(...)=0 helps, for example see the attached file, x = u + v*I

MP_235618.mw

sols = {u[1] = 0., u[2] = 0., u[3] = 0., u[4] = 0., u[5] = 0.,
v[1] = -1.40967543044154,
v[2] = 0.372048922420580,
v[3] = -1.40967543044154,
v[4] = 0.628572434471089,
v[5] = 0.499261752811266}

## reflect the function at the diagonal...

One can "always" plot the inverse function without computing it - just as the reflection at the diagonal:

 > # https://www.mapleprimes.com/questions/235504-Cannot-Plot-Erf-Function
 > plot(erf, -3 .. 3): p1:=%: # remember the plot
 > # reflect the plot at the diagonal, i.e. the line through (x,y) = (0,0) and (1,1) plottools[reflect](p1, [[0,0], [1,1]]); p2:=%:
 >

## caused by "erf"...

I think the problem comes through "erf" for large inputs. You may try the following at the beginning of your sheet Negativity_(v13_beta_gamma_mu):

local erf; #erf(0.1);
cut:=5.8;
erf:= t -> piecewise(t<-cut,-1, t<= cut, :-erf(t), cut <t, 1);
#erf(0.1), erf(30.1), erf(-40.1);

## integral = -1/2*t^(3/2)*Pi^(1/2)*polylog...

You need a formula, which is not implemented (though Maple knows that it is true ...)

Maple Worksheet - Error

Failed to load the worksheet /maplenet/convert/MP_235246.mws .

## well ......

Look up the help, do not use the option "view"

## another presentation...

 > # https://www.mapleprimes.com/questions/235156-Can-We-Find-The-Maximum-Value-Of-This
 > restart;
 > expr:=sin(sqrt(3)*t)*cos(sqrt(3)*t)*(sqrt(3)*cos(sqrt(3)*t) - sin(sqrt(3)*t))/3;
 (1)
 > t_max := '((-arctan(sqrt(RootOf(_Z^3 - 16*_Z^2 + 16*_Z - 3, index = 2))) + Pi))/sqrt(3)'; evalf(%);
 (2)
 > '0=diff(expr,t)': 'eval(%, t =t_max)'; normal(%): convert(%, radical): evala(%);
 (3)
 > # for a compact presentation use an "alias" R = RootOf(_Z^3 - 16*_Z^2 + 16*_Z - 3, index = 2); alias(%): 't_max'=t_max; 'eval(expr, t =t_max)': '%'= simplify(%);
 (4)
 > # just in case one dislikes the sqrt 't_max = ((-arctan(RootOf(_Z^6 - 16*_Z^4 + 16*_Z^2 - 3, index = 2)) + Pi))/sqrt(3)'; convert(%, radical): evala(%): is(%); # true
 (5)
 >

I only post the result since my sheet is too cruddy.

## plot ......

Plotting the rhs lets one guess a singularity in 0 and MultiSeries-series lets me substitute A = B^9 (or even better: A=1/B^9). Then fsolve works (I used simplify assuming positive before that)

It is always a good idea to plot and from that one sees that you run into "numerical singularities". I used afa:=0.277 for the following

plot(X, 0 .. Pi/2); #shows a problem at x towards Pi/2 (using your X modified for T2)

tmp:=eval(sin(x)*T2, x=Pi/2-0.001):    #example

plot(log10(tmp), y=-Pi/6+afa..Pi/6);   #indicates a reason for that

plot(log10(tmp), y=-Pi/128 .. Pi/128); # locates it even more
plot(log10(tmp), y=0.016 .. 0.016005);

For the x-y plane you can guess the trouble for x ~ Pi/2 and y ~ 0 using

'sin(x)*T2';
plot3d(%, x=0 .. Pi/2, y=-Pi/6+afa ..Pi/6, axes=boxed);

'sin(x)*T2'; log10('%');
plot3d(%, x=0 .. Pi/2, y=-0.05 ..0.05, axes=boxed);

'sin(x)*T2'; log10('%');
plot3d(%, x=Pi/2 - 0.05 .. Pi/2.0, y=-0.05 ..0.05, axes=boxed);

I do not have a good idea for that, so you may switch to the rough estimate shown by vv in his contribution. But I wonder how Mathematica treats that peaks and how reliable the result (you forgot to mention it) may be

## "Integral transformations" knows it...

It is a bug I would say. But the package for integral transformations knows it

 > restart: #interface(version);
 > J:=Int(sin(x)*exp(-2*I*Pi*f*x)/x, x = -infinity .. infinity); # upper case Pi
 (1)
 > fourier(F(x),x,t); convert(%, Int); eval(%, F= 'x->sin(x)/x'): eval(%, t=2*Pi*f); # to see what we need as input
 (2)
 > # one has to tell Maple what it should take for Heaviside(0) to write it "piecewise" # here Heaviside(0) = 1 gives the desired (!) result # try  Heaviside(0) = 0 as well for fun ... NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1)):
 > #fourier(sin(x)/x,x,t); eval(%, t=2*Pi*f); fourier(sin(x)/x,x,2*Pi*f);   convert(%, piecewise, f); # note: that implicitly assumes f to be real
 (3)
 >

## ~ 3.515*10^18...

I get 3.515*10^18 (4 leading decimals) within ~ 20 sec by interchanging the integration

Digits:=15;

X:=proc(xx)
local f;
f:=unapply(subs(x=xx, sin(x)*T1), y);
evalf(Int(f,-Pi/6..-Pi/6+afa, epsilon=1e-5, method = _d01ajc));
end proc;

evalf(k*Int(X, 0..Pi/2, method = _d01akc, epsilon = 1e-4));

The problem might stem from oscillations towards Pi/2

MP_234881_text_program.mw

## Use upper case Matrix, not the old style...

Use upper case Matrix, not the old style lower case matrix.

## cleaning up the integrals...

One can simplify the integrals over the complex error function into just one, which speeds up.

I did not use acer's error handling / code (I do not quire understand it), but a cut off.

A proper way for the nasty error  function would be through https://en.wikipedia.org/wiki/Faddeeva_function because of the damping by exp(-kappa^2). There is code for it (through the literature on the German Wikipedia site), but not in Maple.

Edit 13.09.22: the cut off = 6 results from roughly estimating an (absolute) integration error to be less than 1e-15.

 > # https://www.mapleprimes.com/questions/234736-How-To-Reduce-Evaluating-Time
 > restart; interface(version); Digits:=15; local gamma: # or better something like gg instead
 (1)
 > X:=Int(exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*cos(kappa*gamma)/(kappa), kappa = Lambda .. infinity); Y:='Int(-exp(-kappa^2*(delta^2 + 1)/2)*cos(kappa*beta)*Im(exp(kappa*gamma*I)*erf( (gamma+kappa*I)/sqrt(2)))/(kappa), kappa = Lambda .. infinity)'; theJ:= 'sqrt(X^2+Y^2)/2*exp(-alpha^2/2)*lambda^2';
 (2)
 > forget(evalf): 'eval(theJ, [alpha=4,delta=1,Lambda=0.001,beta=8,gamma=3])': CodeTools:-Usage(''%''=evalf(%) ); # this is valid for 15 decimals
 memory used=2.61GiB, alloc change=198.67MiB, cpu time=29.86s, real time=28.76s, gc time=2.45s
 (3)
 > subs(infinity=6, theJ): # cut off numJ2:=subsindets(%,'specfunc(anything,Int)', 'q -> Int(op(q), epsilon = 1e-4, method = _d01akc)'):
 > forget(evalf): 'eval(numJ2, [alpha=4,delta=1,Lambda=0.001,beta=8,gamma=3])': CodeTools:-Usage(''%''=evalf(%) );
 memory used=8.97MiB, alloc change=0 bytes, cpu time=109.00ms, real time=125.00ms, gc time=0ns
 (4)
 >

## plots[display]...

For example:

 > p1:=plot(sqrt(x), x=0..2, color=red): p2:=plot(1/x, x=-3..0, color=blue): plots[display]([p1,p2]);
 >

## suggestion...

Essentially you need Int(x^(1/10)*exp(x),x=0..2) = 6.47950080882301

Int(x^(1/10)*exp(x),x);
IntegrationTools:-Parts(%, exp(x));
value(%):
V:=simplify(%) assuming 0<x;

Plotting "shows" it is continous [Re(V), Im(V)]: plot(%, x=-0.1..2.1);

limit(V, x=2, left)-limit(V, x=0, right);
expand(%);
evalf(%);

"proves" the value
-1/10*(-1)^(9/10)*GAMMA(1/10,-2)+2^(1/10)*exp(2)+1/10*Pi*(-1)^(9/10)/sin(1/10*Pi)/GAMMA(9/10)
is correct

Have not tried for your final integral

## correct call...

You want to write proc(x,y,t), else y is local (as you did) and not known outside

 2 3 4 5 6 7 8 Last Page 4 of 93
﻿