Axel Vogt

5936 Reputation

20 Badges

20 years, 250 days
Munich, Bavaria, Germany

MaplePrimes Activity


These are answers submitted by Axel Vogt

You can try the following (the integrand is rational in C):

EQ1 := convert(eq1, rational); # to avoid floating point numbers

lhs(EQ1);
V1:=value(%) assuming 0<C;
V1:=convert(V1, piecewise, C);

After some time you get an answer. Plot it (in C = 0 ... 1) to see if
you have any chance. You do have a good chance. Howver 'Solve fails
(a bug?). But fsolve will do.

Dito for V2 with C < 0 (I have not checked symmetry).

For V1 or V2=-1140 I get C = +- 0.100907264887074

I think the limit is 1 by the initial value theorem for Laplace transforms.

The 2nd integral works, while the 1st does not. But you forgot to define parameters which is

not n = 1 but n:= 1 (with double point, and so on)

PS: I do not correct it, since handling sheets in document mode is painful


After all the discusion here is the answer (applying Gradshteyn & Ryzhik (2005), Formula 6.512.1 and the Errata for it, giving a a scaled, normalized 2F1, but use the reciprocal Gamma instead of dividing by Gamma):

it is zero, as Markiyan asserted

  Phi:= z -> piecewise(z::nonposint, Limit(1/GAMMA(zeta), zeta=z), 1/GAMMA(z)):

  setting:=mu=1, a=1500,nu=2,b=r; # for Errata
  # versus mu=2, a=r,nu=1,b=1500 for the Book

  # show that the settings give the task
  Int(BesselJ(mu,a*x)*BesselJ(nu,b*x),x = 0 .. infinity):
  % = subs(setting, %);

  # now apply the formula in G&R
  Int(BesselJ(mu,a*x)*BesselJ(nu,b*x),x = 0 .. infinity) =
    b^nu*a^(-nu-1)* GAMMA((mu+nu+1)/2)/GAMMA(nu+1)/GAMMA((mu-nu+1)/2)*
    hypergeom([(mu+nu+1)/2, (mu-nu+1)/2 ],[nu+1],b^2/a^2);

  ``=b^nu*a^(-nu-1)* GAMMA((mu+nu+1)/2)*'Phi'(nu+1)*'Phi'((mu-nu+1)/2)*
    hypergeom([(mu+nu+1)/2, (mu-nu+1)/2 ],[nu+1],b^2/a^2);

  subs(setting, %); eval(%);
  value(%);

                                  = 0

  # and be sure, that the conditions the formula (Errata) are satisfied
  [b < a, 0< a, 0 < b, -1 < Re(mu+nu)]; # Errata: b < a, not a < b !
  subs(setting, %);
  map(is, %) assuming 0 < r, r < 40;

               [b < a, 0 < a, 0 < b, -1 < Re(mu + nu)]
                       [true, true, true, true]

 

Download MP_task_BesselJ.mws

 

Using Maple 17 I would do:
Sum(1/(k^2-p),k=0..infinity); value(%); discont(%, p);

                                    1/2       1/2
                         Pi cot(Pi p   ) p + p
                    -1/2 ------------------------
                                    3/2
                                   p

                                     2
                              {0, _Z1 }

Which says: as far as p is not the square of an integer
Does not look nice or easy to understand. But does it.
gamma is a numerical constant in Maple, you can not assume on it
[solve(eq,beta)];
select('q -> is(0<=q)', %) assuming 1<g;
                               2     1/2
                             (g  - 1)
                            [-----------]
                                  g

Where I would prefer words instead: 1 < q, so sqrt is real & positive
and deviding by a positive does not change it. I see no need to do
work with software here (just my opinion)
 

I guess: k/6/Pi, k = 0 .. 6 and period Pi? But what is the reference for the task?

Edited: I used "convert(%, sin): combine(%): numer(%): collect(%, sin);" for my plot, then it works without increasing numpoints by plot(%, x=-0.1 .. Pi+0.03,  -1 .. 1);

x^t/y^t = (x/y)^t; 
                              t
                             x          t
                            ---- = (x/y)
                              t
                             y

subs(x=1, y=-1, t=1/2, %);
                                -I = I
Here is a different approach, working for larger numbers (but I do not know
to generalize the constraint). If in the following one sets theNumber:=100
then one gets Kitonum's answer (and the constraint is satisfied). But I do
not know, whether the solution is minimal. However it works and does not
suffer from combinatorial explosions. I show it for 200! instead of 100!.

restart: interface(version);
theNumber:=200;
theFactorial:='theNumber!'; `~~`=evalf(%);
                                             375
                     ~ = 0.788657867364791 10


ifactors(theFactorial)[2]:
H:=sort(%);  # prime number decomposition
aF:=nops(%); # Anzahl Faktoren
                               aF := 46

F:=[seq( op(1,H[i]), i = 1 .. aF)]; # Prim-Faktoren
M:=[seq( op(2,H[i]), i = 1 .. aF)]; # Multiplizitäten
K:=[seq( F[i]^M[i],  i = 1 .. aF)]; # primary Komponenten
#'theFactorial=product(K[i], i = 1 .. aF)'; is(%);

r:='r':
sys:=NULL:
for i from 2 to aF - 1 do
  sys:= sys, (K[i]-1) + r[i]*K[i]=(K[i+1]-1) + r[i+1]*K[i+1];
end do:
 
i:=1:
sys:= sys, r[i]*K[i]=(K[i+1]-1) + r[i+1]*K[i+1]:
sys:={sys}: i:='i':

sol:='isolve(sys, t)'; #: expand(%);
sol:=%:
                        sol := isolve(sys, t)

eval(sys, sol): map(is, %);
                                {true}

eval(sys, sol):
X:=lhs( %[1]);  # that is linear in t

0 < X: solve(%): t in evalf(%);

         t  in  RealRange(Open(-0.707274961818205), infinity)

z:=eval(X, t=0); # = coeftayl(X, t=0, 0)

  z := 55779796302805917999479190314474522489737682702230310402429\
        7818973188035586150018567544726120267456829264718819404726\
        7416941228242282110716071584580903149517863596346964322709\
        3217976455085566203844368409127431658463316731949007653535\
        0330428665588959653463722356562158478360867394606691475414\
        2768977006039823370360752097105410374234232308961090041066\
        19977392256259918212890624

`~~`=evalf(z);
                                             375
                     ~ = 0.557797963028059 10

#z/(99!); evalf(%);
evalf(z / (theFactorial/theNumber));

                           141.454992363641

Now verify it is a solution

modp(z*(z+1), theFactorial);
                                  0

For a more detailled look one can use the following code(s):

#for i from 1 to aF do
#  modp(z, K[i])* modp(z+1, K[i]), 'component'[i] = H[i];
#end do; i:='i':

#for i from 1 to aF do
#  modp(z, K[i]), modp(z+1, K[i]), 'component'[i] = H[i]; #, 'elements' = K[i];
#end do; i:='i':

The last code was the reason for the above - feeding it with Kitonum's
solution z and for theFactorial = 100! and experimenting & pondering.

I always hated such notations in Math, but besides Carl's answer:

limit(delta(x), delta = 1);
                                  x

diff(delta(x), x);
                                delta

 

Better multiply by cosh(C) to avoid troubles.

Then note that Maple wants to find a symbolic solution (over the complex numbers) and does not find it (even using 'allvalues' will not help).

So if you want real solution you better use fsolve.

If you plot that task you will see, that there are 2 solutions:

cosh(C+cosh(C)) - 2*cosh(C);
plot(%, C= -2.5 .. 0.5);

fsolve(f, C=-2.5);
                          -2.38320860635885
fsolve(f, C= 0.5);
                          0.323074102021700

convert(%, Diff) will show what you want to see (UpperCase D)

 

@Dima

http://www.mapleprimes.com/view.aspx?sf=202964_Answer/With_par.mw

To be honest I would prefer that you reduce to a more readable question, I do not want to 'crawl' into a procedure. Would you mind to remove the min & max and state in a form like for the "3 integrals"?

I would not say, that I can answer it (or find the time the next days), but it still may be the case, that your Heaviside still means to integrate over a domain (dim=3 in general, yes?). And it may be a way to proceed from there (instead feedeing to 'cuhre' in total).

For short: it is too difficult for me to 'read' your sheet.

PS: I post that as extra 'answer' only to keep the thread more readable.

 

f: X ---> Y, then cokern(f) = Y/image(f) and Y = image(f) (+) complement, where
(+) = direct sum. Hence cokern ~ complement of the image (assume dim < infinity).

To compute it: let Maple compute a base for image(f) and then ask it to enlarge
it for a base of Y. Removing the first should give it. I think that should be
possible with Maple in concrete cases.

Though I never heared that one is asked that in numerical settings and you
want to use convert/rational first.
I get the solution for the system as k2=k1, k3=-2*k1, k3 = 0.158396539759225,
hardware precision.

To give a sketch: as already said write each as one triple integral (removing
min and max) over B1, B2, B3 and use 'convert/rational', Digits=15.

For the above asserted relation one does not need the values. The above gives
B1 - B2 having no parameters anymore and integrating over a = 0 ... 2*Pi can
be done by Maple symbolically. The result is simplified to zero for the ranges
of w and b. Dito for 2*B2 + B3.

So one needs to compute one of those integrals, I did for the 3rd, using B3.

Integrating w.r.t. 'a' is immediate. Now I do it w.r.t. 'w' (because that does
not have 1/sqrt(variable) it seems to be more pleasant). That results in a very
complicated expression (and sometimes Maple seemed to hang up, I lost patients).

Investigating that in detail allows to reduce the range and to eliminate all
the nasty things stemming from the Heaviside function (and to find some
symmetries):

solw7 := -32080347/23129762200*((-160000/117*cos(b)^5+160000/117*cos(b)^3 *
  arcsin(117/200*1/sin(2*b))-10000/117*Pi*(cos(b)-1/2*cos(3*b)-1/2*cos(5*b)) +
  cos(b)*(-20000*cos(4*b)+6311)^(1/2))*Pi/(-49*cos(2*b)+51)^(1/2)

Finally that reduces to compute

4*Int('solw7',
  b=1/2*Pi-arcsin(1/40*634^(1/2)+1/40*166^(1/2)) ..
    1/2*Pi-arcsin(1/40*634^(1/2)-1/40*166^(1/2)), method=_d01ajc );

evalf(%);

                          0.158396539759225

From that one gets k1 and k2.

First 19 20 21 22 23 24 25 Last Page 21 of 93