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

Maple's GE (and GJE) sometimes switches around the rows, but not the columns, so it's unclear what you mean by "...the columns and see where they are after..."

The row switches can be seen by the pivot Vector/Matrix, which is one of the outputs of LUDecomposition (LU is GE in mild disguise). This can be "applied" after the fact to a Vector of names.

acer

Here is a short example which might help you interpret and see the relationship between the algebraic and Matrix forms accepted by LPSolve.

There is no exposed conversion utility (likely because one of the central purposes of the Matrix form is to never form the large expressions in the algebraic form at all). But you might experiment with the internal routine which I use below, in order to see what happens when there are no strict equality constraints, or no inequality constraints, or no simple variable bounds, and so on.

> restart:

> problem := -x-2*y, {7 >= 3*x+5*y, 113 >= -11*x-13*y, 17*x+19*y=23}, x=-29..31, y=-37..41:

> Optimization:-LPSolve(problem);

     [-2.92857142857143, [x = -0.642857142857142, y = 1.78571428571429]]

> kernelopts(opaquemodules=false):

> problem_matrix_form := Optimization:-Convert:-AlgebraicForm:-LPToMatrix(problem)[2..4];

        problem_matrix_form := [-1., -2.], 

          [[-11.  -13.]                               ]                            
          [[          ], [113., 7.], [17.  19.], [23.]], [[-29., -37.], [31., 41.]]
          [[  3.    5.]                               ]                            

> Optimization:-LPSolve(problem_matrix_form);

                  [                   [-0.642857142857142]]
                  [-2.92857142857143, [                  ]]
                  [                   [  1.78571428571429]]

[edited] I might as well just do it. The first two places represent variables ur and dr, and the remaining N places represent w[i], i=1..N.

restart:
randomize():

with(Optimization):
with(Statistics):
with(LinearAlgebra):

N := 500:
R := RandomMatrix(N, outputoptions = [datatype = float[8]]):
ER := Vector[column]([seq(ExpectedValue(Column(R, j)), j = 1 .. N)]):

# algebraic form
st:=time():
W := Vector(N, symbol = w):
S := Vector(N, fill = 1, datatype = float[8]):
z := Multiply(R, Matrix(W)):
con1 := Transpose(W).S = 1:
con4 := seq(z[i][1] - dr >= 0, i = 1 .. N):
con5 := expand(Transpose(W).ER) - ur >= 0:
sol1:=LPSolve(ur+dr, {con1, con4, con5}, seq(w[i]=0..1,i=1..N), maximize = true):
sol1[1];
                        2.68825855741790
time()-st;
                             6.296

# Matrix form
st:=time():
A:=Matrix(N+1,N+2,datatype=float[8]):
A[2..-1,3..-1]:=R:
A[2..-1,1]:=Vector(N,fill=-1,datatype=float[8]):
A[1,3..-1]:=ER:
A[1,2]:=-1:
MatrixScalarMultiply(A,-1,inplace):
b:=Vector(N+1,datatype = float[8]):
#A,b;
Aeq:=Matrix(1,N+2,fill=1,datatype=float[8]):
Aeq[1,1..2]:=Vector(2,datatype=float[8]):
beq:=Vector(1,[1],datatype=float[8]):
#Aeq,beq;
c:=Vector(N+2,[1,1],datatype=float[8]):
#c;
bl:=Vector(N+2,datatype=float[8]):
bl[1..2]:=Vector(2,fill=-infinity):
bu:=Vector(N+2,fill=1,datatype=float[8]):
bu[1..2]:=Vector(2,fill=infinity,datatype=float[8]):
#bl,bu;
sol2:=Optimization:-LPSolve(c,[A,b,Aeq,beq],[bl,bu],maximize = true):
sol2[1];
                        2.68825855741790
time()-st;
                             3.735

So, half the time, for N=500. There are also memory savings on the order of the size of all the constraints as symbolic expressions, which I didn't bother to measure. There might be a relatively small bit more time saved, using ArrayTools routines like BlockCopy and Fill.

acer

Obviously we cannot be completely sure about the cause of the problem, if we cannot see the full worksheet that exhibits the problem.

But I suspect that Georgios K. is on the right tack, when looking at the incomplete worksheet containing just the final formulae.

I suspect roundoff as the underlying cause because, as Georigios mentions, the plots stabilize when computed at a higher working precision. I can actually get `function1` to plot in what seems to be a reliable and accurate manner with Digits set only to 33, provided that I first manipulate `function1` with the `evalc` command. That command will treat `theta` as purely real and under that assumption will simplify the formula enough to remove all merely ostensible occurences of `I` the imaginary unit. Not only does evalc(function1) appear to plot more accurately without such a high Digits setting as `function1` requires, but it also does so measurably faster.

The other reason I suspect roundoff is because the order of terms in PROD and SUM dags still depends on memory value for ordering. As I tried to illustrate in this old blog post, the order in which multiplicands are stored internally in a product data structure can greatly affect floating-point evaluation. The same is true of sums (which I hope to blog about soon, though I must apologize, for I have written that before...). The upshot is that session dependent ordering of arithmetic terms can magnify and enhance numerical difficulties at a given working precision. How high to increase Digits, the working precision, in order to overcome, is yet another tough topic.

Again, without the Original Poster's full worksheets, we cannot investigate and discover the underlying cause in a definitive way. But some of the familiar bells that are ringing seem to include roundoff and memory-based ordering of terms.

To trlewis, please feel free to contact me if you would consider my investigating your worksheets in a confidential way.

Since simplification (including using evalc, say) need not necessarily improve the numeric robustness of any formula to be evaluated in floating-point, then for now I'd simply concur with Georgios: increase Digits.

Example_alt01.mw Roundoff may be more severe in your original full worksheets than in this modified sheet, depending on term ordering issues. It seems like a possibility.

I do not know what, if anything, MuPAD might try to do, automatically, to accomodate numerical issues for such examples.

acer