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

parse(sprintf("%.2f",654321.3678918));

                              654321.37

I can not remember who initially posted that trick.

I can not confirm this, neither in older Maple 15 nor in concurrent version 17
(Win 32, Maple 32 classical interface), may you post details and upload a sheet

http://www.mapleprimes.com/ViewTemp.ashx?f=492_1388940326/zero_vector_matri.mws

completesquare(%, k0)

I would do, but do not like to re-type code posted as images

I guess the following happens: the integrand is 'sent as is' to IEEE evaluation.

But that can not recognize the damping and evaluates each of the factors, which
finally avoids results 'good' enough.

So - as a rule - one has to prepare the task for numerics. Always.

In this case it is easy: just combine into one exp, that is

  Digits:= trunc(evalhf(Digits)):
  uu := x -> Int(exp(-(xi-x)^2+xi), xi = -infinity .. infinity, method = _d01amc);
  evalf(uu(0.));

               2.27587579446875

That is even fast enough to plot the function.


PS: one better uses option epsilon first (digits only reduces precision)
in case one needs such.

PPS: you can check by using a symbolic solution,

  Int(exp(-(xi-x)^2+xi), xi = -infinity .. infinity); value(%);

                                         1/2
                          exp(x + 1/4) Pi

I avoid it to modify inputs (and if your x is assigned to a value, then it
will fail to work).

May be the following old post from the usenet helps a bit (have not tried,
whether Google provides the complete thread):

--------------------------------------------
Subject: Re: passing by reference
Date: 15 Jun 2003 21:49:43 GMT
From: israel@math.ubc.ca (Robert Israel)
Organization: ITServices, University of British Columbia
Newsgroups: comp.soft-sys.math.maple

In article <h05Ha.206462$ro6.5581866@news2.calgary.shaw.ca>,
Elaine Jackson <elainejackson7355@home.com> wrote:
>How do you pass variables by reference to a Maple procedure? TIA.

Maple's parameter passing is normally "call by evaluated name".
Thus if you have a procedure f and a variable v whose value is y,
if you say f(v) then the actual parameter passed to f is y.  To pass
the name v, there are basically two methods:

1) Delay evaluation of the parameter using quotes, i.e. f('v').
2) Define f so its parameter will not be evaluated, i.e.
 f:= proc(p::uneval) ...

Robert Israel                                israel@math.ubc.ca
Department of Mathematics        http://www.math.ubc.ca/~israel

--------------------------------------------

You remember n*(n+1)/2 to count such? Try it with n = 1, 2, 3 to see it.

Then n*(n+1)/2=2013; fsolve(%); tells you to take n = 63 and now check
(look up the help for the operator "$" and try it with 2$2; 3$3; etc).


[seq(n$n, n=1 .. 63)]; 
%[2013];
nops(%%);

If f and g are inverse to each other you only know that their compositions are the identity.

"It" does not have problems, but the usage is not what you want:

  convert(ExcelTools:-Import("sample.xls", "Sheet1", "B2:B10"),list)

works correctly, giving a list of 9 floats. However then writing

  S:= {op(...)};

makes it a set (same as convert/set). Which 1) eliminates duplicates
and 2) destroys the ordering (seems it gets numerical sorted).

Have not tried to understand Carl's original code.

The first 2 are quadratic equations, the 3 other ones are trigonomtric equations.

I played with it, solving the trigonometrics first (hoping Maple finds enough, then checking your assumptions) and then feeding the other. Without checking everything

What is the solution you already know?

int( 1/(x+sqrt(x)),x); convert(%, ln); simplify(%);
limit(%, x=1, left);
limit(%%, x=1, right);

# both are correct

2*ln(1+sqrt(x));; diff(%, x); 
1/(x+sqrt(x)) = simplify(%); is(%);

ln(x-1)+2*arctanh(x^(1/2)); diff(%, x); 
1/(x+sqrt(x)) = simplify(%); is(%);

# see the difference

2*ln(1+sqrt(x)); subs(x=x+y*I, %); plot3d(abs(%), x=-1..1, y=-1..1, axes=boxed);

ln(x-1)+2*arctanh(x^(1/2)); plot3d(abs(%), x=-1..1, y=-1..1, axes=boxed);

First: double precision (in the sense of IEEE 754 ---> Wikipedia) always means 32 Bit, whatever you have as 64 on your machine as OS and SW

Second: You are comparing 3 SW systems, which probably handle that differently. Plus it will depend on the setting for the system and the commands.

Third: What you see may not what you have. You see decimals, while double prec is binary and you see a display given by a conversion routine. The same for the inputs, the must be 'translated' to binary (=base 2) from decimals (=base 10).

Set a:=916.536, b:=3.3. What is displayed as a+b ?

Matlab probably will convert a to some IEEE value A, same for B, the will add both and finally displays it as decimal.

MMA probably will convert, add and display as well. But here it will depend on the settings for working and display precision (and perhaps on some automatics which only an experienced user may know - not me).

Maple will do similar, but here I am aware of many things, so I can say a bit about mixing base 2 = evalhf and base 10:

Be aware, that you want Digits:=18 for displaying/reading double precision and decimal representation.

Maple will give a+b as 919.836. Applying evalhf tells Maple to convert that to a IEEE value, which is 4045481510587859/4503599627370496 * pow(2,10) and then computing & displaying it in decimals, giving 919.836000000000012.

For Matlab the command will be different: it probably *first* convert a and b to   8061935914202431/9007199254740992 * pow(2,10),   3715469692580659/4503599627370496 * pow(2,2) then adds them - howver not add exactly, but according to double precision giving

  8090963021175717/9007199254740992 * pow(2,10)

and from that a (no, not unique) decimal representation is

  919.835999999999898

MP_200300.mws is the worksheet

plot( x -> Int(exp(-(x-xi)^2)*f(xi),xi = -infinity .. infinity, method = _d01amc), -10 .. 10);