Axel Vogt

5821 Reputation

20 Badges

20 years, 221 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

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

 

Occasionally I do the following

f := x -> piecewise(0 < x,-1,x <= 0,1): # double point = quiet
'f(x)': '%'=%;                          # a simple way to display such things
                             { -1        0 < x
                      f(x) = {
                             { 1         x <= 0

 

One can read it as a(i+1) = i^2 - a(i), which is a recursion.
Maple can solve this, i.e. the term a(k) is given "directly":

  restart;
  rsolve(a(i)+a(i+1)=i^2, a(k)); # see the help for the command

                        k
               a(0) (-1)  - 1 + (k + 1) (k/2 + 1) - 2 k

 


For specific values one has to provide a(0), of course.

In "standard interface" one can insert a table (through the menu bar) and each cell can have execution groups to display the plots.

If you dislike to show the commands the you can hide them (via menu/show) - but (for me) that only works for the entire worksheet.

The correct way ist

x^(5/3)-5*x^(2/3);
[Re(%), Im(%)];
plot(%, x=-3..7, color=[red, blue]);

Instead of using floating point numbers, like 0.3, one can use decimal fraction, like 3/10.

Then the DE can be solved symbolically. Plotting indicates a zero at r ~ 1.5*10^5  1.5*10^6 and that complex values might occure.

Analyzing and using my own zero finder (like bracketting) gives a zero between 1337252.92663784 and 1337252.92663785


Be aware that this depends on the exactness of your input parameters (i.e. how many decimal places are used) and how many Digits are used for computations.

I used Digits:=15

 

Int(1/u*t*exp(-t/u), t = 0 .. infinity);
value(%) assuming 0<u;

Numerically one splits in Pi/2 * k, k integer, because of needed 'smoothness' ~ subdivisions. One would need k ~ 10^6 for double precision, I think.

Manual support gives Sum(2*Si(Pi*(2*k+1))-Si(2*k*Pi)-Si((2*k+2)*Pi)+(Si(k*Pi)-Si((k+1)*Pi))*(-1)^k,k = 0 .. infinity) and that numerically gives the expected result.

First 8 9 10 11 12 13 14 Last Page 10 of 92