Remove the two extra spaces. One space is between DEplot and the ( bracket, and the other space is between diff and the ( bracket.

Those spaces make Maple think that those are (implicit) multiplications DEplot*(...) and diff*(...) when you enter it in 2D Math input mode (the default entry mode).

acer

The LSSolve result is not necessarily higher. It's higher because you chose a value of ER=10, which is used only to fill y the target for the least-squares fit. It fits to target data, by minimizing error. Give it an ad hoc target, and get an ad hoc result. Why shouldn't the LSSolve result depend on your arbitrary ER value? 

I'm not suggesting that you obtain data2 in order to figure out what ER ought to be, as a general methodology. (That would ve absurd, as data2 is an answer to your problem.)

restart:
interface(warnlevel=0):
randomize(): 
with(ListTools): 
with(LinearAlgebra): 
with(ArrayTools): 
with(Statistics): 
with(plots): 
with(Optimization): 

nstock := 400: 
nr := 50: 
##ER := 10: # arbitrary target will lead to arbitrary result. 10 -> too high on avg.

W := Vector([seq(w[i], i = 1 .. nstock)]): 
R := RandomMatrix(nr, nstock, outputoptions = [datatype = float[8]]): 
Cov := CovarianceMatrix(R): 
ev := Vector([seq(ExpectedValue(Column(R, i)), i = 1 .. nstock)], datatype = float[8]):

s2 := Optimization[QPSolve](Transpose(W).Cov.W-Transpose(W).ev):
data2 := eval(R.W, s2[2]):

magic:=Variance(data2);

                   1.2392001406827267

ER := 2*magic;
                   2.4784002813654533

y := Vector(nr, fill = ER, datatype = float[8]): 

s1 := Optimization[LSSolve](convert(R.W-y, list)):
data1 := eval(R.W, s1[2]):

display({CumulativeSumChart(data1, color = red, legend = "LSSolve",
                            markers = false, thickness = 3),
         CumulativeSumChart(data2, color = green, legend = "QPSolve",
                            markers = false, thickness = 3)});

acer

You could try and run the loop in Maple itself, using the SetSubstitutions and Simulate commands? The runtime of MapleSim is generally accessible in the Maple interface, from a set of commands (an API).

Basic syntax might look something like this:

A := MapleSim:-LinkModel('filename' = "yoursavedmodel.msim"):

A:-GetSubstitutions(); # shows you the present set of parameter values

for x from value1 to value2 by increment_value do

   A:-SetSubstitutions({`modelandparamname` = newvaluecomputedfromx})

   arrayofresults := A:-Simulate(output=datapoint, anyotheroptions);

   # You may want returnTimeData=true as an option to Simulate().
   # Each column of output Array is from a separate probe (by row w.r.t. time)

   ExportMatrix(filename, Matrix(arrayofresults));

   # Or export with writedata, or whatever.
   # You might want to use a different export file for each value of x,
   # by using something like cat("basename","_",convert(x,string)) for 'filename'.

end do;

I suspect that instead of calling SetSubstitutions() you might be able to pass the changing parameter value directly to Simulate() as an option like 'params={name=value}'

acer

You should be able to do this in at least two ways.

One way is to not declare intXdA as an Array, so that assigning to its entries implictly make it a table. In Maple, tables don't have a fixed number of positions. Then it gets converted to an Array after the computation of all its entries.

Another way is to create the Array only after n is known.

restart:

flsp:= proc(f,ug,og)
local i, n, y, intXdA;
   y:=unapply(f,x):
   if 0>ug then
      if og=0 then n:=1
      elif og>0 then n:=2;
      end if;
   end if: 
   for i from 1 by 1 to n do
      # intXdA is a table, at this moment
      intXdA[i]:=evalf(int(x*y(x),x=ug..og)):
   end do;
   # convert the table into an Array
   Array(1..n,intXdA);
end:

flsp(x,-4,2);
                                 [24., 24.]

flsp(x,-4,0);
                                [21.33333333]

restart:

flsp:= proc(f,ug,og)
local i, n, y, intXdA;
   y:=unapply(f,x):
   if 0>ug then
      if og=0 then n:=1
      elif og>0 then n:=2;
      end if;
   end if:
   # Now that we know n, create the Array
   intXdA:=Array(1..n):
   for i from 1 by 1 to n do
      intXdA[i]:=evalf(int(x*y(x),x=ug..og)):
   end do;
   intXdA;
end:

flsp(x,-4,2);
                                 [24., 24.]

flsp(x,-4,0);
                                [21.33333333]

I made intXdA the return value of the procedure, which is a nicer way to program than to make it a global. Note the use of square-brackets [] when indexing into either table or Array.

You wrote that you have Maple 7. I didn't check this in Maple 7. In modern Maple an Array is "growable" if instead using round-brackets () to access entries, and that allow Array intXdA to be created only as size 1..1 ie. even before n is known.

acer

Most or all of the locals will be assigned values of either integers or floats. If you pass in argument Lint as an Array to hold the integers, and argument Lfloat an Array to hold the floats, then you might be able to get rid of most or all locals and replace them with references to entries of those two Arrays.

Sure, that makes you have to manage the extra bookkeeping. (I've never measured for any performance impact.) Apologies if you've already considered this idea.

acer

Since fsolve will only try a limited number of initial starting points, finding one whch converges can sometimes depend upon supplying ranges for the variables which are not "too" wide.

In the example, it is easy enough to find ranges which work for all seven entries of S (ie. values of w), by looking at the solutions for the few values of w for which convegence succeeds with the Poster's original ranges.

restart:

wvals:=[1.2,1.6,2,2.5,3,3.5,4]:
S:=[seq({w=wvals[i]},i=1..7)];

          [{w = 1.2}, {w = 1.6}, {w = 2}, {w = 2.5}, {w = 3}, {w = 3.5}, 

           {w = 4}]

sys := {exp(-.1204819277*(2.4039*t+15.44745000*t^2-11.03552334*t^3+2.595300000*t^4+.508258/t-44.6834-(2.40397*ln(t)+30.8949*t-16.55325000*t^2+3.460399999*t^3+.2541195000/t^2-49.812126)*t)/t)*a2*a4-a1*a3, exp(-.1204819277*(-2.071844454/t+.3999293136/t^2+2.897999999*t^3-0.2404762368e-1/t^3-6.278629824*t^2+1.49670934*t-.7274250000*t^4+4.532401680*ln(1000*t)+4.532401680*ln(298)+134.6934679-(-4.532276160/t-1.035931727/t^2-.9699000000*t^3+.2666044800/t^3+4.347000000*t^2-12.55719488*t-0.1794936000e-1/t^4-27.04489066*ln(1000*t)+27.04489066*ln(298)+28.54167*ln(t)+190.6774129)*t)/t)*a4*((1/2)*a1+(1/2)*a2+(1/2)*a3+(1/2)*a4+(1/2)*a5)-a1*a2, (1/4)*exp(-.1204819277*(95.3768*t-71.65195000*t^2+24.69369999*t^3-3.579500000*t^4+1.105339/t+179.76736-(95.37701*ln(t)-143.3039*t+37.04055000*t^2-4.772666666*t^3+.5525695000/t^2+363.377422)*t)/t)*(a1+a2+a3+a4+a5)^2*a5*a4-a1^3*a2, 2*a1+2*a4+4*a5-1.6-2*w, a2+a3+a5+a6-1, a2+2*a3+a4-.77-w, -202.86-180.476*w-a1*(33.0661*t-5.681700000*t^2+3.810933333*t^3-.6932000000*t^4+.158558/t-9.9807)-a2*(25.5675*t+3.048050000*t^2+1.351533333*t^3-.6678250000*t^4-.1310/t-118.0118)-a3*(24.9973*t+27.59345000*t^2-11.23045667*t^3+1.987075000*t^4+.1366/t-403.5951)-a4*(30.092*t+3.416250000*t^2+2.264466667*t^3-.6336000000*t^4-0.821e-1/t-250.8806)-a5*(-.703*t+54.23865000*t^2-14.17383333*t^3+1.465675000*t^4-.678565/t-76.84066)-a6*(27.04489066*t+.2287298242*t^2-4.532401680*ln(1000*t)+2.181502454/t-.3999293136/t^2+0.2404762368e-1/t^3-11.80536790-4.532401680*ln(298))}:

st:=time():
sol:='sol':
for i from 1 to nops(S) do
  sol[i]:=fsolve(eval(sys,S[i]),{a1=0..2,a2=0..2,a3=0..2,a4=0..10,a5=0..1,a6=0..2,t=0..2});
  if not type(eval(sol[i],1),specfunc(anything,fsolve)) then
     # i, max.abs. error, solution
     print([i, max(map(abs,evalf(eval(eval(sys,S[i]),sol[i])))), sol[i]]);
  end if;
end do:
time()-st;

 [         -11                                          
 [1, 5.1 10   , {a1 = 0.4745463021, a2 = 0.05174213138, 

   a3 = 0.1964321579, a4 = 1.525393553, a5 = 0.00003007250511, 

                                      ]
   a6 = 0.7517956383, t = 1.049462130}]
 [         -7                                          
 [2, 1.2 10  , {a1 = 0.5998025943, a2 = 0.07298483347, 

   a3 = 0.2484467941, a4 = 1.800121578, a5 = 0.00003791366446, 

                                      ]
   a6 = 0.6785304588, t = 1.060375417}]
 [         -7                                          
 [3, 1.5 10  , {a1 = 0.7251005506, a2 = 0.09503840192, 

   a3 = 0.3000768144, a4 = 2.074807969, a5 = 0.00004574008213, 

                                      ]
   a6 = 0.6048390436, t = 1.067956587}]
 [       -8                                         
 [4, 7 10  , {a1 = 0.8817218880, a2 = 0.1233029268, 

   a3 = 0.3642649879, a4 = 2.418167097, a5 = 0.00005550730265, 

                                      ]
   a6 = 0.5123765780, t = 1.074696743}]
 [       -8                                        
 [5, 4 10  , {a1 = 1.038324469, a2 = 0.1520588490, 

   a3 = 0.4281980719, a4 = 2.761545007, a5 = 0.00006526187271, 

                                      ]
   a6 = 0.4196778172, t = 1.079572257}]
 [          -7                                        
 [6, 1.29 10  , {a1 = 1.194906156, a2 = 0.1811351045, 

   a3 = 0.4919605334, a4 = 3.104943829, a5 = 0.00007500746271, 

                                      ]
   a6 = 0.3268293547, t = 1.083262425}]
 [         -8                                        
 [7, 1.3 10  , {a1 = 1.351469280, a2 = 0.2104315203, 

   a3 = 0.5556036263, a4 = 3.448361227, a5 = 0.00008474652456, 

                                      ]
   a6 = 0.2338801069, t = 1.086152488}]
                             4.150

One can also run this at higher working precision,

Digits:=20:

st:=time():
sol:='sol':
for i from 1 to nops(S) do
  sol[i]:=fsolve(eval(sys,S[i]),{a1=0..2,a2=0..2,a3=0..2,a4=0..10,a5=0..1,a6=0..2,t=0..2});
  if not type(eval(sol[i],1),specfunc(anything,fsolve)) then
     # i, max.abs. error, solution
     print([i, max(map(abs,evalf(eval(eval(sys,S[i]),sol[i])))), sol[i]]);
  end if;
end do:
time()-st;

[         -17                                
[1, 1.3 10   , {a1 = 0.47454630227099735486, 

  a2 = 0.051742131398933455701, a3 = 0.19643215794129892897, 

  a4 = 1.5253935527184686864, a5 = 0.000030072505266979386860, 

                                                         ]
  a6 = 0.75179563815450063595, t = 1.0494621299751955770}]
[         -17                                
[2, 2.6 10   , {a1 = 0.59980259451150934909, 

  a2 = 0.072984833487173778850, a3 = 0.24844679417683001007, 

  a4 = 1.8001215781591662010, a5 = 0.000037913664662224945514, 

                                                         ]
  a6 = 0.67853045867133398614, t = 1.0603754166160084808}]
[         -17                                
[3, 1.1 10   , {a1 = 0.72510055082673461650, 

  a2 = 0.095038401944926777977, a3 = 0.30007681452327544115, 

  a4 = 2.0748079690085223397, a5 = 0.000045740082371521882767, 

                                                         ]
  a6 = 0.60483904344942625899, t = 1.0679565867975334120}]
[         -17                                
[4, 1.3 10   , {a1 = 0.88172188825734010461, 

  a2 = 0.12330292680996198531, a3 = 0.36426498802662595140, 

  a4 = 2.4181670971367861119, a5 = 0.000055507302936891753254, 

                                                         ]
  a6 = 0.51237657786047517153, t = 1.0746967427846299290}]
[         -17                               
[5, 1.6 10   , {a1 = 1.0383244694955892758, 

  a2 = 0.15205884905244901302, a3 = 0.42819807209461689471, 

  a4 = 2.7615450067583171976, a5 = 0.000065261873046763296416, 

                                                         ]
  a6 = 0.41967781697988732897, t = 1.0795722568764767296}]
[          -17                               
[6, 1.86 10   , {a1 = 1.1949061567313915910, 

  a2 = 0.18113510449275594600, a3 = 0.49196053358241408129, 

  a4 = 3.1049438283424158914, a5 = 0.000075007463096258800517, 

                                                         ]
  a6 = 0.32682935446173371391, t = 1.0832624249502111602}]
[          -17                               
[7, 2.64 10   , {a1 = 1.3514692802731791716, 

  a2 = 0.21043152031278464376, a3 = 0.55560362650519318177, 

  a4 = 3.4483612266768289927, a5 = 0.000084746524995917846612, 

                                                         ]
  a6 = 0.23388010665702625662, t = 1.0861524880905071629}]
                             4.430

With a guess as to the nature of the problem (ranges too wide for a limited number of initial starting points) it is easier to find out that only the original range for a5 need be changed in order to find solutions for all seven w values.

Digits:=10:

st:=time():
sol:='sol':
for i from 1 to nops(S) do
  sol[i]:=fsolve(eval(sys,S[i]),{a1=0..10,a2=0..10,a3=0..10,a4=0..10,a5=0..0.1,a6=0..10,t=0..10});
  if not type(eval(sol[i],1),specfunc(anything,fsolve)) then
     # i, max.abs. error, solution
     print([i, max(map(abs,evalf(eval(eval(sys,S[i]),sol[i])))), sol[i]]);
  end if;
end do:
time()-st;

 [         -11                                          
 [1, 5.1 10   , {a1 = 0.4745463021, a2 = 0.05174213138, 

   a3 = 0.1964321579, a4 = 1.525393553, a5 = 0.00003007250511, 

                                      ]
   a6 = 0.7517956383, t = 1.049462130}]
 [         -7                                          
 [2, 1.2 10  , {a1 = 0.5998025943, a2 = 0.07298483347, 

   a3 = 0.2484467941, a4 = 1.800121578, a5 = 0.00003791366446, 

                                      ]
   a6 = 0.6785304588, t = 1.060375417}]
 [         -7                                          
 [3, 1.5 10  , {a1 = 0.7251005506, a2 = 0.09503840192, 

   a3 = 0.3000768144, a4 = 2.074807969, a5 = 0.00004574008213, 

                                      ]
   a6 = 0.6048390436, t = 1.067956587}]
 [       -8                                         
 [4, 7 10  , {a1 = 0.8817218880, a2 = 0.1233029268, 

   a3 = 0.3642649879, a4 = 2.418167097, a5 = 0.00005550730265, 

                                      ]
   a6 = 0.5123765780, t = 1.074696743}]
 [       -8                                        
 [5, 4 10  , {a1 = 1.038324469, a2 = 0.1520588490, 

   a3 = 0.4281980719, a4 = 2.761545007, a5 = 0.00006526187271, 

                                      ]
   a6 = 0.4196778172, t = 1.079572257}]
 [          -7                                        
 [6, 1.29 10  , {a1 = 1.194906156, a2 = 0.1811351045, 

   a3 = 0.4919605334, a4 = 3.104943829, a5 = 0.00007500746271, 

                                      ]
   a6 = 0.3268293547, t = 1.083262425}]
 [         -8                                        
 [7, 1.3 10  , {a1 = 1.351469280, a2 = 0.2104315203, 

   a3 = 0.5556036263, a4 = 3.448361227, a5 = 0.00008474652456, 

                                      ]
   a6 = 0.2338801069, t = 1.086152488}]
                             9.843

To be sure, fsolve could benefit from an additional option to specify the total number of loop iterations and also the total number of starting points, as well options for tolerances and working precision.

acer

> alpha:=0.618:

> sprintf( "The value of alpha is %a.", alpha );
                 "The value of alpha is .618."

acer

The control character of a tab can be printed by using an escaped character with backslash (\).

fprintf("c:/temp/text.txt","%a\t%a",A,B):

fclose("c:/temp/text.txt");

And now that file text.txt contains,

A	B

with a tab between them.

See ?backslash and ?printf

acer

h:=(1+i)^2+(1+i)^3*A+(1+i)^4*B;

                      2          3            4  
               (1 + i)  + (1 + i)  A + (1 + i)  B

algsubs(1+i=g,h);

                         2      3      4
                        g  + A g  + B g 

Of course, in this simple example subs will also work, but one has to "solve" for `i` (mentally works, in this example),

subs(i=g-1,h);

                         2      3      4
                        g  + A g  + B g 

This kind of attempt with subs is also possible here (automating the "solving"),

subs(isolate(i+1=g,i),h);

                         2      3      4
                        g  + A g  + B g 

acer

> restart:

> h:=Statistics:-Distribution(PDF=(x->Dirac(x-b))):

> eval(h[':-PDF']);

                               x -> Dirac(x - b)

> map(FromInert,indets(remove(type,[op(ToInert(eval(h[':-PDF'])))],
>             specfunc(anything,{_Inert_OPTIONSEQ,_Inert_LOCALSEQ,
>                                _Inert_PARAMSEQ})),
>               specfunc(string,_Inert_NAME)));

                                  {Dirac, b}

You should be able to sieve that result, or the inert body, to remove names of applied function calls (like `Dirac`).

acer

It looks like Maple 9 is confused about what to do with the structure of the float object itself.

> kernelopts(version);

            Maple 9.03, IBM INTEL LINUX, Oct 1 2003 Build ID 141050

> dismantle([10000.]); # writedata works ok for this list         

LIST(2)
   EXPSEQ(2)
      FLOAT(3): 10000.
         INTPOS(2): 10000
         INTPOS(2): 0

> dismantle(map(x->1/x, [.0001])); # writedata fails for this list   

LIST(2)
   EXPSEQ(2)
      FLOAT(3): .1e5
         INTPOS(2): 1
         INTPOS(2): 4

These both seem to do "better", in 9.03,

> restart:

> B:=[ 3.5, 0.0001, 73.45]; A:=[map(x->1.0*1/x,B)];

                           B := [3.5, 0.0001, 73.45]

               A := [[0.2857142857, 10000.00000, 0.01361470388]]

>  writedata("a.txt",A); # ok      
                
>  writedata("a.txt",A,float); # ok

and produce,

% more a.txt
0.2857142857	10000	0.01361470388

It's interesting. The printf and fprintf commands also seem to malfunction here.

> B:=[ 3.5, 0.0001, 73.45]; A:=[[seq(1/op(j,B), j=1..3)]];

                           B := [3.5, 0.0001, 73.45]

                 A := [[0.2857142857, 10000., 0.01361470388]]

> printf("%{n}.11g\n", Matrix(A));       
                 
0.2857142857 1 0.01361470388

>  B:=[ 3.5, 0.0001, 73.45]; A:=[map(x->1.0*1/x,B)];    
  
                           B := [3.5, 0.0001, 73.45]

               A := [[0.2857142857, 10000.00000, 0.01361470388]]

> printf("%{n}.11g\n", Matrix(A));          
        
0.2857142857 10000 0.01361470388

It looks like this was fixed in the next major release, 9.5.

acer

There are several ways to determine whether two lines are perpendicular. Which is easiest to use can depend on which particular form one has for the lines. See here as the "Graph of functions" section for two easy ways, both of which should require no deep thinking to do in Maple.

acer

Is this just a QP minimization problem in disguise, rather than an NLP minimization problem? By which I mean, a quadratic objective and linear constraints?

Sure, your con3 constraint is nonlinear, involving a 2-norm. But it is the only place that variable `r` appears in constraints, and your objective is simply `r`. So can't you just make the norm part of con3 as a new objective, and then get rid of con3 as a constraint?

Note that the setup of this problem gets expensive quickly, as `nstock` grows. And there is a lot of room for improvements to the first parts which generate the randomized data. So, unless all that is just a stand-in for some other means of data acquisition then you can get a bit of speedup just by refining all that code. I mean, everything up to the Cholesky decomposition of `Cov`. On the machine used for the timings below, at nstock=200 the construction of `Cov`, `R` and the constraints takes about 70 seconds. What comes before constructing `Cov` can likely be reduced to under 4-5 seconds.

I set nstock=200 and compared these two Optimization calls below. I'll leave out all the set up here, except for the constraint assignments,

con1 := add(W[i], i = 1 .. nstock) <= 1:
con2 := EV.W-pr >= 0:

NRW2 := Norm(R.W, 2, conjugate = false):
con3 := NRW2 - r <= 0:

# As you solved it with NLPSolve...

sol1 := CodeTools:-Usage(NLPSolve(r,{con1,con2,con3},seq(w[i]=0..1,i=1..nstock)));

memory used=1.74GiB, alloc change=55.49MiB, cpu time=60.68s, real time=60.70s

sort(select(t -> rhs(t)>0 and lhs(t)<>r, sol1[2]));

        [w[21] = 0.0270237593173089, w[25] = 0.0322056751304918, 

          w[35] = 0.000847852100225448, w[67] = 0.0492087037211113, 

          w[78] = 0.0438242963135371, w[86] = 0.0416518751919245, 

          w[88] = 0.0417036500369194, w[89] = 0.0717631323305688, 

          w[100] = 0.00837739752003451, w[107] = 0.137179654344253, 

          w[109] = 0.107669160523004, w[125] = 0.0289124908386128, 

          w[131] = 0.160793105594057, w[153] = 0.0236924836513717, 

          w[157] = 0.0341339386566807, w[180] = 0.0119978162730577, 

          w[185] = 0.0405857945037750, w[190] = 0.129027000241401, 

          w[192] = 0.00940221371168578]

r = eval(r, sol1[2]);

                            r = 2.03629907765324

# And now...

sol2 := CodeTools:-Usage(QPSolve(NRW2^2,{con1,con2},seq(w[i]=0..1,i=1..nstock)));

memory used=448.67MiB, alloc change=0 bytes, cpu time=7.71s, real time=7.72s

sort(t -> rhs(t)>0 and lhs(t) <> r, sol2[2]));

        [w[21] = 0.0270237589337289, w[25] = 0.0322056750319796, 

          w[35] = 0.000847852090717340, w[67] = 0.0492087042316293, 

          w[78] = 0.0438242958978819, w[86] = 0.0416518740369907, 

          w[88] = 0.0417036488838085, w[89] = 0.0717631314373649, 

          w[100] = 0.00837739703763107, w[107] = 0.137179652525067, 

          w[109] = 0.107669163079282, w[125] = 0.0289124902420816, 

          w[131] = 0.160793108650386, w[153] = 0.0236924831211481, 

          w[157] = 0.0341339380256316, w[180] = 0.0119978167550648, 

          w[185] = 0.0405857946727649, w[190] = 0.129027001732921, 

          w[192] = 0.00940221361388309]

r = sqrt(eval(NRW2^2, sol2[2])); # `r` in the original sense of `con3`

                            r = 2.03629907765382

Unless I did something quite wrong, those two attempts produce essentially the same solution, but the QPSolve is 7 times faster.

And the QPSolve use can likely also be improved using Matrix form for its input.

acer

> evalc(ln(x)) assuming x<0;

                         ln(-x) + I Pi

acer

Thy this code by Joe Riel.

acer