Axel Vogt

5018 Reputation

19 Badges

15 years, 315 days
Munich, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

"If the file is being opened for writing and does not exist (and is not already open), it is created, if possible; otherwise, fopen generates an error.  If the file does exist (and is not already open), it is truncated. Writing starts at the beginning of the file if WRITE was specified, and at the end of the existing file if APPEND was specified."

Then you can use iostatus(), for example.

In my ancient times of coding at least a check for existence was done (and available handles), though I do not expect Maple to provide all the stuff for a coding laguage.


 

Digits:=15:
f:= x ->
  evalf( Int('t -> Re(exp(LambertW(1/(-1+t))*(-1+t)))', 1 .. x, method=_d01ajc) )
 +evalf( Int('t -> Im(exp(LambertW(1/(-1+t))*(-1+t)))', 1 .. x, method=_d01ajc) )*I:

g:= (-1)/f+1;

[Re@g, Im@g]:
plot(%, 0 .. 2, -10 .. 10, color=[red,blue, grey]);

-1/f+1

 

 

 


 

Download 222149.mw

 

It is a convention what 0^0 should mean and Maple treats it as follows:

`0^0`: %= parse(%);

                               0^0 = 1
`0^r`: %= parse(%);

                               0^r = 0
`0^1`: %= parse(%);

                               0^1 = 0
`0^(-1)`: %= parse(%);

   Error, numeric exception: division by zero

Actually you expect your result to be 11 * `0^0` = 11.


You can however use your own powering:

pow:= (t,n) -> piecewise( t=0 and n=0, 1, t^n);

m:=2; # or your m=10
x:=0;

Sum(f(x,j)*(m-j+1), j=0..m); # or your original double sum
value(%);
subs(f=pow, %);
eval(%);

                          2
                        -----
                         \
                          )   f(0, j) (3 - j)
                         /
                        -----
                        j = 0


                   3 f(0, 0) + 2 f(0, 1) + f(0, 2)


                3 pow(0, 0) + 2 pow(0, 1) + pow(0, 2)


                                  3

 

Int(sin(phi)*sin(-arctan((p*cos(phi)+Z)/sqrt((p*sin(phi))^2+r^2)))*hypergeom([5/2, 3/2], [5/2], sin(-arctan((p*cos(phi)+Z)/sqrt((p*sin(phi))^2+r^2)))^2)/((p*sin(phi))^2+r)^2, p = -1 .. 1);
tmp:=subs(phi=Pi/6, %); # = your function, evaluated
plot(op(tmp)); # "visually" correct

 

simplify(tmp);
   -32*3^(1/2)*(Int(p*(p^2+25)/((p^2+100)^(3/2)*(p^2+20)^2), p = -1 .. 1))

Now just argue by symmetry in zero.

 

 

 

Use fsolve( ..., complex)

As far as I can see you have a set of 36 quadratic equations. It is very likely that there is more than 1 solution (and it is not clear at all why those [complex ?] solutions are discrete). It is also possible that there is NO solution (avoiding "evalf" seems to indicate that).

PS: you can avoid the lengthy display by just uploading you sheet or a pdf of it.

rsolve({a(n) = 2*a(n-1) + 4*a(n-2), a(0)=1, a(1)=2}, a(k));

Your integrand has a singularity in 0.161144379880543e-3, so you can not use numerical integration directly. However the integral exists as Cauchy Principal Value (if it is what you want):

(E*sigma*(Int((-xi^2+1)*(-KImax^2+Kfc^2)/((KImax^2-Ksc^2)*((-R^4+1)*(KImax^2+Ksc^2)+eta*E*sigma)),
  l = l0 .. b, CPV))/alpha0);

value(%);

    8491.79571174594

where b is your upper bound 12.1*10^(-3)

For me it is a bit difficult to understand your task and quite personal notations (and you upload is not the same if I execute your sheet, try to use Maple sheets with a restart at the beginning).

May be you want the following

tmp6:=(1/2)*erf((1/2)*sqrt(2)*d6)+1/2;

'Int(exp(-beta^2/T)* tmp6 /sqrt(Pi*T), beta = 0 .. infinity)';

value(%) assuming 0 < T:
simplify(%);
                     1    /1  (1/2)   \   1
                     - erf|- 2      d6| + -
                     4    \2          /   4

 

I convert to rational and 'optimize' the computation - usually Maple should care for rounding issues (checking in details may be difficult) - and get 0.29531564 quickly.

If you are interested in more than just this specific value it may be worth to go such a route (else it would not matter that much if a programm needs 1 sec or 1 h).

Perhaps you re-state your (original) task so the structure can be seen - which can rarely done having the code alone.

MP221823.mws

MP221823.pdf

 

simplify(ii_inf)  assuming 1/2<s,s<1;
MultiSeries:-asympt(%,x,1) assuming 1/2<s,s<1;


  1/8*4^s*Pi^(1/2)*GAMMA(s-1/2)/sin(Pi*s)/GAMMA(s)+O((1/x)^(2*s-1))

I would proceed similar, but a bit different (except you only care for the 'length' - but why?).

Usually one expects such a fraction to be 'exact' for the precision one needs.

Thus I would do

v := 1433162862156581221897903/4340366287316864960103999;
oldDigits:= Digits:
  Digits:=10:        # fall down to desired precision
  evalf(v):
  vNew:=convert(%, rational);
Digits:=oldDigits:   # restore settings


                         1433162862156581221897903
                    v := -------------------------
                         4340366287316864960103999


                                    15251
                            vNew := -----
                                    46188

In this example vNew = v is correct up to at most the last decimal place
in a precision of 10 decimal places.

# symmetric in x and then integrand has an inverse
ann:= 0<x,0<q,1<y;

Int(sech(x)^q, x);
IntegrationTools:-Change(%, x=arccosh(y), y);
combine(%) assuming ann;
simplify(%) assuming ann;
value(%) assuming ann;
A:=eval(%, y = cosh(x)) assuming ann;

 

# result:

-1/q*cosh(x)^(-q)*hypergeom([1/2, 1/2*q],[1+1/2*q],1/(cosh(x)^2));

Your f is abs(g) and g is positive over the range 1 ... 9. So just write Int(g, p= 1 .. 9), which numerically evaluates quickly (and even faster if you use 'simplify' before feeding it to evalf).


 

#VectorCalculus:-int(sin(x^2)*cos(y^2), [x, y] = Circle(`<,>`(0, 0), 1), inert);

sin(x^2*cos(y)^2)*cos(x^2*sin(y)^2)*x; #int(%, y);
4*Int(%, y=0..Pi/2);
IntegrationTools:-Change(%, y = arcsin(t), t); value(%);
Int(%, x=0..1);
value(%);
evalf(%);

sin(x^2*cos(y)^2)*cos(x^2*sin(y)^2)*x

 

4*(Int(sin(x^2*cos(y)^2)*cos(x^2*sin(y)^2)*x, y = 0 .. (1/2)*Pi))

 

-4*x*(Int((sin(x^2*t^2)*cos(-x^2)+cos(x^2*t^2)*sin(-x^2))*cos(x^2*t^2)/(-t^2+1)^(1/2), t = 0 .. 1))

 

x*sin(x^2)*Pi

 

Int(x*sin(x^2)*Pi, x = 0 .. 1)

 

(1/2)*Pi-(1/2)*cos(1)*Pi

 

.722091449378415

(1)

 


 

Download MP-220842.mw

Shift it to 0, find an ODE for it, solve it and determine parameters:

task:=hypergeom([1/3, 2/3],[3/2],27/4*z^2*(1-z));
h:=eval(%, z=1+x);

# find differential equation
f(x)=h;
PDEtools[dpolyform](%, no_Fn):
ode:=op(1, %)[1];
dsolve(ode):
H:=rhs(%);

# check ode
odetest(f(x)=H,ode);
odetest(f(x)=h,ode);

# find parameters
'eval(h, x=1)'; convert(%, StandardFunctions): evala(%):
c1:=%=eval(H, x=1); #convert(%, StandardFunctions): evala(%);

'eval(h, x=2)'; convert(%, StandardFunctions): evala(%):
c2:=%=eval(H, x=2);

param:=solve({c1,c2});

# re-substitute
eval(H, param);
eval(%, x=z-1);

 

2 3 4 5 6 7 8 Last Page 4 of 86