This seems almost as fast, when increasing to mm=7 and the range of the dsolving to 0.2 shows it to be slower (64bit Linux, 15.01).

And then for mm=8 and higher this seems to take much longer (I didn't wait, after 3-4 times the duration as the other).

Earlier code, also using evalf[10] for the numerical integration all done beforehand,

restart:
st:=time():
# Maple15\xsys11demo.mw, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure
with(StringTools):
FormatTime(%c);
with(DEtools):with(plots):
printlevel := 1:
mm := 7: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

for n to nn do
for nd to nn do
for md to mm-1 do
integ[n,nd,md]:=evalf[10](Int( sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)
                           *sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2),
                           y=0..a,x=0..a))
end do
end do
end do;

dproc3 := proc (nn, t, Y, YP) local n, nd, md;
  for n to nn do
    YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]
             *add(add(Y[nd]*conjugate(Y[nd])*integ[n,nd,md],
                      md = 1 .. mm-1), nd = 1 .. nn);
  end do;
end proc:

ic3 := array([seq(1/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:

p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0.,
            initial = ic3, procvars = dvars3, range=0..0.2, method=rkf45):

FormatTime(%c);
odeplot(p,[seq([t,abs(c[n](t))],n=1 .. nn)]): # suppressed display only
time()-st;

                   "Fri Mar 30 14:57:50 2012"
                   "Fri Mar 30 14:57:52 2012"
                             1.770

Newer code, for same parameters. Sometimes it takes 1.86 sec on my machine, and sometimes 3.3 sec.

restart:
st:=time():
# Maple15\xsys11demo.mw, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure
with(StringTools):
FormatTime(%c);
with(DEtools):with(plots):
printlevel := 1:
mm := 7: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

dproc3 := proc (nn, t, Y, YP)
local n, nd, md;
   for n to nn do
      YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]*add(add(Y[nd]*conjugate(Y[nd])
         *evalf[10](Int(sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)
            *sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2),
         [y=0..a,x=0..a])), md=1..mm-1), nd=1..nn);
   end do;
end proc:

ic3 := array([seq(1/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:

p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0.,
            initial = ic3, procvars = dvars3, range=0..0.2, method=rkf45):

FormatTime(%c);
odeplot(p,[seq([t,abs(c[n](t))],n=1 .. nn)]): # suppressed display only
time()-st;

                   "Fri Mar 30 14:58:54 2012"
                   "Fri Mar 30 14:58:56 2012"
                             1.860

Changing to mm=8 and the latter code above seems to "go away".

Did I copy or use it not as intended?

acer

This seems almost as fast, when increasing to mm=7 and the range of the dsolving to 0.2 shows it to be slower (64bit Linux, 15.01).

And then for mm=8 and higher this seems to take much longer (I didn't wait, after 3-4 times the duration as the other).

Earlier code, also using evalf[10] for the numerical integration all done beforehand,

restart:
st:=time():
# Maple15\xsys11demo.mw, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure
with(StringTools):
FormatTime(%c);
with(DEtools):with(plots):
printlevel := 1:
mm := 7: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

for n to nn do
for nd to nn do
for md to mm-1 do
integ[n,nd,md]:=evalf[10](Int( sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)
                           *sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2),
                           y=0..a,x=0..a))
end do
end do
end do;

dproc3 := proc (nn, t, Y, YP) local n, nd, md;
  for n to nn do
    YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]
             *add(add(Y[nd]*conjugate(Y[nd])*integ[n,nd,md],
                      md = 1 .. mm-1), nd = 1 .. nn);
  end do;
end proc:

ic3 := array([seq(1/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:

p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0.,
            initial = ic3, procvars = dvars3, range=0..0.2, method=rkf45):

FormatTime(%c);
odeplot(p,[seq([t,abs(c[n](t))],n=1 .. nn)]): # suppressed display only
time()-st;

                   "Fri Mar 30 14:57:50 2012"
                   "Fri Mar 30 14:57:52 2012"
                             1.770

Newer code, for same parameters. Sometimes it takes 1.86 sec on my machine, and sometimes 3.3 sec.

restart:
st:=time():
# Maple15\xsys11demo.mw, standard worksheet, BoltzQM, Shift & enter, ADD not SUM, 1-D math for procedure
with(StringTools):
FormatTime(%c);
with(DEtools):with(plots):
printlevel := 1:
mm := 7: nn := 3: a := 5.: e := 1.:  R := .1: h := 1.:

dproc3 := proc (nn, t, Y, YP)
local n, nd, md;
   for n to nn do
      YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]*add(add(Y[nd]*conjugate(Y[nd])
         *evalf[10](Int(sin(Pi*n*x/a)*sin(Pi*nd*x/a)*sin(Pi*n*y/a)
            *sin(Pi*nd*y/a)/sqrt((x-y)^2+4*R^2*sin(Pi*md/mm)^2),
         [y=0..a,x=0..a])), md=1..mm-1), nd=1..nn);
   end do;
end proc:

ic3 := array([seq(1/sqrt(nn), n = 1 .. nn)]):
dvars3 := [seq(c[n](t), n = 1 .. nn)]:

p := dsolve(type=numeric, number = nn, procedure = dproc3, start = 0.,
            initial = ic3, procvars = dvars3, range=0..0.2, method=rkf45):

FormatTime(%c);
odeplot(p,[seq([t,abs(c[n](t))],n=1 .. nn)]): # suppressed display only
time()-st;

                   "Fri Mar 30 14:58:54 2012"
                   "Fri Mar 30 14:58:56 2012"
                             1.860

Changing to mm=8 and the latter code above seems to "go away".

Did I copy or use it not as intended?

acer

@casperyc You have changed the methodology, and should beware of premature evaluation.

Try this,

restart:

myf := Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=-infinity..0, epsilon=1e-5, method=_d01amc)
     + Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=0..infinity, epsilon=1e-5, method=_d01amc):

newmyf:=student[changevar](x=u*sigma,myf,u):

evalf(eval(myf,[mu=-1.645074,sigma=15000])),
evalf(eval(newmyf,[mu=-1.645074,sigma=15000])); # they agree, good

myf0:=unapply('evalf'(myf),[mu,sigma]):
myf0(-1.645074,16500); # does not succeed

newmyf0:=unapply('evalf'(newmyf),[mu,sigma]):
newmyf0(-1.645074,16500); # but this does succeed

plot('newmyf0'(-1.645074,sigma),sigma=17000..18000); # without quotes, you will wait!

plot('newmyf0'(-1.645074,sigma),sigma=1000..18000);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=100..18000,axes=box);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=3000..7000,axes=box);

@casperyc You have changed the methodology, and should beware of premature evaluation.

Try this,

restart:

myf := Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=-infinity..0, epsilon=1e-5, method=_d01amc)
     + Int(exp(-(x-mu)^2/2/sigma^2)/sigma/sqrt(2*Pi)/(1+exp(-x)),
           x=0..infinity, epsilon=1e-5, method=_d01amc):

newmyf:=student[changevar](x=u*sigma,myf,u):

evalf(eval(myf,[mu=-1.645074,sigma=15000])),
evalf(eval(newmyf,[mu=-1.645074,sigma=15000])); # they agree, good

myf0:=unapply('evalf'(myf),[mu,sigma]):
myf0(-1.645074,16500); # does not succeed

newmyf0:=unapply('evalf'(newmyf),[mu,sigma]):
newmyf0(-1.645074,16500); # but this does succeed

plot('newmyf0'(-1.645074,sigma),sigma=17000..18000); # without quotes, you will wait!

plot('newmyf0'(-1.645074,sigma),sigma=1000..18000);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=100..18000,axes=box);

plot3d('newmyf0'(mu,sigma),mu=-10..10,sigma=3000..7000,axes=box);

@Jarekkk One of the benefits of having Mat_proc_newer accept Matrix Mf and Vector yf as float[8] rtable arguments is that it can be called under evalhf. That's what makes Minimize work fast with it.

And it can also be used that way for single calls, eg,

restart:
NN:=7054:
M := Matrix(NN, 60):
for i to 10 do M[1, i] := 1 end do:
y:=[0.8045, 0.1834, 0.0006]:

Mat_proc_newer := proc(x,M,NN,y)
  local i::integer, j::integer;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:

Mf:=Matrix(M,datatype=float[8]):
yf:=Vector(y,datatype=float[8]):

CodeTools:-Usage(Mat_proc_newer(0.3,Mf,NN,yf));
memory used=79.00MiB, alloc change=11.75MiB, cpu time=1.50s, real time=1.50s
                      0.4419552141658895

CodeTools:-Usage(evalhf(Mat_proc_newer(0.3,Mf,NN,yf)));
memory used=0.58KiB, alloc change=0 bytes, cpu time=265.00ms, real time=263.00ms
                      0.441955214165889509

@Jarekkk One of the benefits of having Mat_proc_newer accept Matrix Mf and Vector yf as float[8] rtable arguments is that it can be called under evalhf. That's what makes Minimize work fast with it.

And it can also be used that way for single calls, eg,

restart:
NN:=7054:
M := Matrix(NN, 60):
for i to 10 do M[1, i] := 1 end do:
y:=[0.8045, 0.1834, 0.0006]:

Mat_proc_newer := proc(x,M,NN,y)
  local i::integer, j::integer;
  for i from 2 to NN do
    M[i, 1] := M[i-1, 2];
    M[i, 60] := M[i-1, 59];
    for j from 2 to 59 do
      M[i, j] := M[i-1, j] + x * (M[i-1, j-1] - 2*M[i-1, j] + M[i-1, j+1])
    end do
  end do;
return add(((1/10)*add(M[NN, i+10*j], i = 1 .. 10)-y[j])^2, j = 1 .. 3);
end proc:

Mf:=Matrix(M,datatype=float[8]):
yf:=Vector(y,datatype=float[8]):

CodeTools:-Usage(Mat_proc_newer(0.3,Mf,NN,yf));
memory used=79.00MiB, alloc change=11.75MiB, cpu time=1.50s, real time=1.50s
                      0.4419552141658895

CodeTools:-Usage(evalhf(Mat_proc_newer(0.3,Mf,NN,yf)));
memory used=0.58KiB, alloc change=0 bytes, cpu time=265.00ms, real time=263.00ms
                      0.441955214165889509

@herclau The original submission has Qnum as an expression, and not a procedure or operator. (I had forgotten that I'd earlier accomodated this using `unapply`, in my answer above.)

So, using Maple 15.01 (64bit, Windows 7), a result better than that from the discrete plot grid is obtained just using,

> Optimization:-Maximize(Qnum, Do = 0.5e-1 .. .5, N = 1 .. 800);

        [24560.2479846227870, [Do = 0.32224998719071896,  N = 142.62810464042423]]

@herclau The original submission has Qnum as an expression, and not a procedure or operator. (I had forgotten that I'd earlier accomodated this using `unapply`, in my answer above.)

So, using Maple 15.01 (64bit, Windows 7), a result better than that from the discrete plot grid is obtained just using,

> Optimization:-Maximize(Qnum, Do = 0.5e-1 .. .5, N = 1 .. 800);

        [24560.2479846227870, [Do = 0.32224998719071896,  N = 142.62810464042423]]

@Thomas Michael Curson We cannot see the embedded image in your post. Is it an image of more code, perhaps pasted (from a Mac?) instead of uploaded using the green up-arrow in the posting editor?

@maple fan It turns out that `int` can also be used here, since it will try numerical intergation if the range values (x below) are floats (or if something like fsolve wraps evalf calls around it, when evaluating at specific points). With Maple 15.01 (Linux, 64bit),

> st:=time():

> fsolve(x->int('g(t)', t=0..x)-5, 0..10);

                                  7.500000000

> time()-st;

                                     0.390

@maple fan It turns out that `int` can also be used here, since it will try numerical intergation if the range values (x below) are floats (or if something like fsolve wraps evalf calls around it, when evaluating at specific points). With Maple 15.01 (Linux, 64bit),

> st:=time():

> fsolve(x->int('g(t)', t=0..x)-5, 0..10);

                                  7.500000000

> time()-st;

                                     0.390

And also see Preben Alsholm's procedure, ResizePlot, here.

And also see Preben Alsholm's procedure, ResizePlot, here.

@hermitian As you've guessed, CodeTools:-Usage is just something that you can wrap around a computational function call, to get an easy printing of the time and memory resources that get used. It's not necessary here, and might not have been available in Maple 13.

The `epsilon` is an option to evalf(Int(...)) the numeric integration (quadrature) command. It is the accuracy tolerance. You see, fsolve is picky and might not return a result unless it is satisified with the convergence of its result. And fsolve might try and attain an accuracy that depends upon the current Digits setting at the level at which it was called. So I set the quadrature accuracy-tolerance `epsilon` to be small enough to get 10 digits of accuracy, but put that inside the EnergyDensity procedure which had Digits raised high enough (hopefully) to allow evalf/Int to satisfy such a tolerance.

The other stuff inside evalf(Int(...)) was to force use of a method that might be suitable for a semi-infinite range of numeric integration. See the evalf/Int help-page. If Digits=15 had not been high enough to allow evalf/Int to get the provided `epsilon` accuracy then I'd have had to set Digits even higher within `EnergyDensity` and not specified _d01amc which only works up to Digits=15 and hardware double precision.

acer 

@hermitian As you've guessed, CodeTools:-Usage is just something that you can wrap around a computational function call, to get an easy printing of the time and memory resources that get used. It's not necessary here, and might not have been available in Maple 13.

The `epsilon` is an option to evalf(Int(...)) the numeric integration (quadrature) command. It is the accuracy tolerance. You see, fsolve is picky and might not return a result unless it is satisified with the convergence of its result. And fsolve might try and attain an accuracy that depends upon the current Digits setting at the level at which it was called. So I set the quadrature accuracy-tolerance `epsilon` to be small enough to get 10 digits of accuracy, but put that inside the EnergyDensity procedure which had Digits raised high enough (hopefully) to allow evalf/Int to satisfy such a tolerance.

The other stuff inside evalf(Int(...)) was to force use of a method that might be suitable for a semi-infinite range of numeric integration. See the evalf/Int help-page. If Digits=15 had not been high enough to allow evalf/Int to get the provided `epsilon` accuracy then I'd have had to set Digits even higher within `EnergyDensity` and not specified _d01amc which only works up to Digits=15 and hardware double precision.

acer