Hopefully I'm now understanding your explanation of the problem.

If you are trying to find N and non-zero X such that K.X=0 hold approximately, then Matrix K should have an eigenvalue close to zero. So, you could try and find an N such that the smallest absolute value of K's eigenvalues is close to zero.

restart:
randomize():
with(LinearAlgebra):
Digits:=20:
interface(warnlevel=0):

# eps is a tolerance, how close to zero for an accepted eigenvalue
eps := 1e-10:

m := 3: # increase, while adding other "known" Matrices
A:=LinearAlgebra:-RandomMatrix(m,generator=-1.0..1.0):
B:=LinearAlgebra:-RandomMatrix(m,generator=-1.0..1.0):
C:=LinearAlgebra:-RandomMatrix(m,generator=-1.0..1.0):

# routine `z` finds N such that K(N) has near-zero as an eigenvalue
z:=N->min(abs(LinearAlgebra:-Eigenvalues(A+B*N+C*N^2))):

K:=N->A+B*N+C*N^2:

sol:=Optimization:-Minimize(z,-600..700):

success:=false:
if sol[1] > eps then
   error "no value of N found giving a near-zero eigenvalue";
else
   success:=true;
   Nsol := sol[2][1];
   Ksol := K(Nsol);
   (evals,Xtry) :=LinearAlgebra:-Eigenvectors(Ksol);
   Xsol:=Matrix([seq(`if`(eps>=abs(evals[i]),Xtry[1..-1,i],NULL),i=1..m)]);
end if:

if success then
  Nsol, Xsol;
  K(Nsol);
  K(Nsol) . Xsol;
end if;

                                   [0.27797677482440519410 + 0. I ]
                                   [                              ]
          -0.48910186887266487015, [-0.69924889690086809645 + 0. I]
                                   [                              ]
                                   [-0.65861968755962728970 + 0. I]

 [ 0.06944242311757242511   0.95596328018536999674  -0.98562630390977603059]
 [                                                                         ]
 [-0.07922533330327931080  -0.32370935082657659106   0.31024065587712641287]
 [                                                                         ]
 [ -1.0251409778044306607  -0.44163425638354885242  0.036207365656109929698]

                         [              -12       ]
                         [-1.53402370 10    + 0. I]
                         [                        ]
                         [             -12        ]
                         [3.85882681 10    + 0. I ]
                         [                        ]
                         [              -12       ]
                         [3.634613295 10    + 0. I]

If you want extra qualifications on N (eg. that N>0 say) then change the range option in the call to `Minimize`, such as 0..infinity or 0..100 etc.

We don't know what your known Matrices are. The above will not succeed for every random A, B, and C;

acer

If I understand your question then you could put `option remember` on the routines which compute the numerical solution. See this comment (and those that follow it).

acer

It is harder to find because it is a Comment. And a Comment is not an Anwser, Question, or Post. And there is no easy way to see a list of your Comments, from your Profile page.

But here is a link to all your postings.

acer

It is easy enough to rotate around another vertical line (what the Asker calls a z-axis) by shifting the origin in the x-y plane though a simple change of variables in the original expression.

After that, the plotted axes (using the usual axes system, thus allowing for tickmarks, placement, etc) have their origin at the desired new point.

If the new origin is to be the minimum of a paraboloid, then the `Minimize` function can find that directly. Setting up derivatives and solving is not necessary.

Tickmarks which show the shifted origin are easy to contruct.

restart:

P := 10.715691225750079453-.52581617253900429073*x+.50951826882694897598*x^2
     -32.921695149953028411*y+45.846724612790471380*y^2-5.9292825745917680635*x*y:

sol := Optimization:-Minimize(P);

          sol := [-0.7380131370346, [x = 4.17647809408785609, y = 0.629109433919071370]]

xsol,ysol := op([2,1,2],sol), op([2,2,2],sol);

                   xsol, ysol := 4.17647809408785609, 0.629109433919071370

xwidth,ywidth := 20, 5:
maxwidth := max(xwidth,ywidth):

newP := subs(x=x-xsol,y=y-ysol,P): # shift origin, in formula P

plotP := plot3d(newP,x=-xwidth..xwidth,y=-ywidth..ywidth):

numframes:=25:

plots:-animate(plottools:-rotate
               , [plotP,A*2*Pi/numframes,[[0,0,0],[0,0,1]]]
               , A=1..numframes
               , axes=normal, labels=[x,y,z]
               , tickmarks=[[seq(i-xsol=i,i=trunc(-maxwidth-xsol)
                                                  ..trunc(maxwidth+xsol),
                                 maxwidth/5)],
                            [seq(i-ysol=i,i=trunc(-maxwidth-ysol)
                                                  ..trunc(maxwidth+ysol),
                                 maxwidth/5)],
                            default]);

This also runs OK in Maple 13.02, Standard GUI.

acer

You show lists, but you mention Matrices.

For lists, you can concatenate using Kitonum's way.

For Matrices, it could be entered like so,

A := <<1|2>>;

                                 A := [1  2]

B := <<3|4>>;

                                 B := [3  4]

C := <A|B>;

                              C := [1  2  3  4]

acer

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