A (less flashy) alternative is,

a:=12+10*x^4+7*x^3-89*x^2+65*x:

r:=[fsolve(a)];

     -3.624808208, -0.1524813627, 1.095553748, 1.981735823

seq(eval(a,x=t),t=r);

                    -7       -9      -8      -7
                5 10  , -7 10  , 3 10  , 1 10  

A (less flashy) alternative is,

a:=12+10*x^4+7*x^3-89*x^2+65*x:

r:=[fsolve(a)];

     -3.624808208, -0.1524813627, 1.095553748, 1.981735823

seq(eval(a,x=t),t=r);

                    -7       -9      -8      -7
                5 10  , -7 10  , 3 10  , 1 10  

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.