Advanced Search

Axel Vogt

5906 Reputation

20 Badges

19 years, 305 days
Munich, Bavaria, Germany
Contact Axel Vogt

MaplePrimes Activity


These are answers submitted by Axel Vogt

for example...

Axel Vogt 5906
January 26 2023
0 4

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

Axel Vogt 5906
January 07 2023
2 2

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:=%:

 
> 

 

Download MP_35504_plot_inverse_function.mw

caused by "erf"...

Axel Vogt 5906
December 30 2022
0 2

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

Axel Vogt 5906
November 26 2022
4 2

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 .

Download MP_235246.mws

well ......

Axel Vogt 5906
November 19 2022
0 1

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

another presentation...

Axel Vogt 5906
November 13 2022
2 0


 

> 

# 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/3)*sin(3^(1/2)*t)*cos(3^(1/2)*t)*(3^(1/2)*cos(3^(1/2)*t)-sin(3^(1/2)*t))

 (1)
> 

t_max := '((-arctan(sqrt(RootOf(_Z^3 - 16*_Z^2 + 16*_Z - 3, index = 2))) + Pi))/sqrt(3)'; evalf(%);

(-arctan(sqrt(RootOf(_Z^3-16*_Z^2+16*_Z-3, index = 2)))+Pi)/sqrt(3)

  

1.39084587051934

 (2)
> 

'0=diff(expr,t)': 'eval(%, t =t_max)';
normal(%):
convert(%, radical):
evala(%);

eval(0 = diff(expr, t), t = t_max)

  

0 = 0

 (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(%);

R = RootOf(_Z^3-16*_Z^2+16*_Z-3, index = 2)

  

t_max = (1/3)*(-arctan(R^(1/2))+Pi)*3^(1/2)

  

eval(expr, t = t_max) = (1/3)*(3^(1/2)*R^(1/2)+R)/(1+R)^(3/2)

 (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

t_max = (-arctan(RootOf(_Z^6-16*_Z^4+16*_Z^2-3, index = 2))+Pi)/sqrt(3)

  

true

 (5)
> 

 

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

Download MP_235156.mw

 

plot ......

Axel Vogt 5906
October 12 2022
0 3

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)

T2 / another answer...

Axel Vogt 5906
October 05 2022
0 9

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

Axel Vogt 5906
October 01 2022
3 2

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

> 

# https://www.mapleprimes.com/questions/234913-Int-Gives-Wrong-Answer-For-Improper-Integral
> 

restart: #interface(version);
> 

with(inttrans):
#? inttrans # help about this package
> 

J:=Int(sin(x)*exp(-2*I*Pi*f*x)/x, x = -infinity .. infinity); # upper case Pi

Int(sin(x)*exp(-(2*I)*Pi*f*x)/x, x = -infinity .. infinity)

 (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

fourier(F(x), x, t)

  

Int(F(x)*exp(-I*x*t), x = -infinity .. infinity)

  

Int(sin(x)*exp(-(2*I)*Pi*f*x)/x, x = -infinity .. infinity)

 (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

Pi*(-Heaviside(f-1/(2*Pi))+Heaviside(f+1/(2*Pi)))

  

piecewise(f < -(1/2)/Pi, 0, f < (1/2)/Pi, Pi, (1/2)/Pi <= f, 0)

 (3)
> 

 

Download MP_234913.mw

~ 3.515*10^18...

Axel Vogt 5906
September 25 2022
1 5

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

Axel Vogt 5906
September 23 2022
2 0

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

cleaning up the integrals...

Axel Vogt 5906
September 12 2022
2 4

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

`Standard Worksheet Interface, Maple 2021.2, Windows 7, November 23 2021 Build ID 1576349`

  

15

 (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';

Int(exp(-(1/2)*kappa^2*(delta^2+1))*cos(kappa*beta)*cos(kappa*gamma)/kappa, kappa = Lambda .. infinity)

  

Int(-exp(-(1/2)*kappa^2*(delta^2+1))*cos(kappa*beta)*Im(exp(I*kappa*gamma)*erf((gamma+I*kappa)/sqrt(2)))/kappa, kappa = Lambda .. infinity)

  

(1/2)*sqrt(X^2+Y^2)*exp(-(1/2)*alpha^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

  

eval(theJ, [alpha = 4, delta = 1, Lambda = 0.1e-2, beta = 8, gamma = 3]) = 0.730427393719220e-3*lambda^2

 (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

  

eval(numJ2, [alpha = 4, delta = 1, Lambda = 0.1e-2, beta = 8, gamma = 3]) = 0.730427393783615e-3*lambda^2

 (4)
> 

 

Download MP_234736.mw

plots[display]...

Axel Vogt 5906
August 29 2022
2 0

For example:


 

> 

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

 
> 

 


 

Download MP_234721.mw

suggestion...

Axel Vogt 5906
August 02 2022
2 3

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

Axel Vogt 5906
July 03 2022
2 0

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