What do you expect to get, by asling for the Jacobian of a scalar valued expression such as your obj??

Are you perhaps thinking instead of the Gradient? So your end goal is the Hessian? The Jacobian of the gradient?

Why not use the Hessian command directly, if that's what you're after?

So, after defining obj and loading Student:-VectorCalculus, either,

Hessian(obj,..);

or,

Jacobian(Gradient(obj,..),..);

including any extra arguments wanted.

What kind of result are you hoping for, for this particular example though?

acer

If you are trying to find the general form of the eigenvectors (as columns of your `NewMatrix3`) of the Matrix of data in your `NewInput3` then your methodology looks suspect.

I do not see how your general equation,

{seq(seq((NewMatrix3 . NewInput3(1..-1,i))[j]=(v[i]* NewInput3(1..-1,i))[j], j=1..3), i=1..3)}

properly represents the definition

A . X[i] = lambda[i] . X[i]

where NewInput3 = A is the data, NewMatrix3[..,i] = X[i] is the ith eigenvector, and v[i] = lambda[i] is the ith eigenvalue. Using those equivalences then shouldn't you instead have,

NewInput3 . NewMatrix3[..,i] = v[i] . NewMatrix3[..,i]

?

As for how LinearAlgebra:-Eigenvectors does it for a Matrix A of exact rational data, the rough answer is that it first finds the eigenvalues. And then for each eigenvalue v[i] it evaluates the characteristic Matrix A-lambda*Id at lambda=v[i] and then computes the nullspace (kernel) of that.

For floating-point data the eigen-solving is passed to a LAPACK function. If you have float data then you might want to consider letting Maple use the fact that NewInput3 = InputMatrix3^%T . InputMatrix3 is symmetric. In the case of float data that will allow Maple to produce purely real results with no very small imaginary (numerical roundoff) artefacts. Ie,

NewInput3 := Matrix( InputMatrix3^%T . InputMatrix3, shape=symmetric );

Or you could consider squaring the singular values of InputMatrix3, and taking care for the ordering. Consider,

restart:
with(LinearAlgebra):

InputMatrix3 := Matrix([[31.25,30.8,30.5],[30.8,30.5,0],[30.5,0,0]]):

map(t->t^2, SingularValues(InputMatrix3));

                             [4789.06261226748]
                             [                ]
                             [591.210617040193]
                             [                ]
                             [284.319270692331]

Eigenvalues(Matrix(InputMatrix3^%T . InputMatrix3, shape=symmetric));

                             [284.319270692331]
                             [                ]
                             [591.210617040193]
                             [                ]
                             [4789.06261226748]

Eigenvalues(InputMatrix3^%T . InputMatrix3);

                          [4789.06261226748 + 0. I]
                          [                       ]
                          [591.210617040192 + 0. I]
                          [                       ]
                          [284.319270692331 + 0. I]

So,...

G := Matrix(InputMatrix3^%T . InputMatrix3, shape=symmetric):

evals:=Eigenvalues(G);

                                  [284.319270692331]
                                  [                ]
                         evals := [591.210617040193]
                                  [                ]
                                  [4789.06261226748]

seq( NullSpace(CharacteristicMatrix(G,evals[i])), i=1..3 );

     /[ 0.326650403460829]\    /[-0.588317874678543]\    /[0.739717238039373]\ 
     |[                  ]|    |[                  ]|    |[                 ]| 
    < [-0.737693385300777] >, < [ 0.330570994334320] >, < [0.588669081053437] >
     |[                  ]|    |[                  ]|    |[                 ]| 
     \[ 0.590853605559242]/    \[ 0.737973506325627]/    \[0.326017055932820]/ 

Eigenvectors(G);

   [284.319270692331]  [-0.326650403460830  -0.588317874678544  0.739717238039373]
   [                ]  [                                                         ]
   [591.210617040193], [ 0.737693385300778   0.330570994334321  0.588669081053437]
   [                ]  [                                                         ]
   [4789.06261226748]  [-0.590853605559242   0.737973506325627  0.326017055932820]

For exact data Maple will compute the nullspace by solving the linear system with the zero-Vector as the right hand side. For floats it will use... (wait for it..), the full SVD.

Why do you want to do it all by hand?

acer

Did you intend a separator in front of "Users", so that it was fully qualified as a file location?

Ie, "/Users/..."

If you don't, then it'll try to use a location relative to whatever the Maple command `currentdir()` returns. So, with your original syntax it might possibly have succeeded and produced a file with some odd long name.

acer

The distinction is mostly between interative applications written with Maplets and those written with Embedded Components.

The Math Apps are examples of the latter, as are the Mobius Apps. The last few releases of Maple have included a growing number of such Math Apps, as featured items (see here and here).

Embedded components reside in a .mw Worksheet or Document. These kinds of apps can also be run interactively in the free Maple Player, as well as in the Maple Standard GUI. And these kinds of worksheet apps done with embedded components have been placed on both the regular and the Mobius Cloud. The Explore command produces a kind of inlined worksheet interative application done with embedded components, and this recent Appclication Center item contains notes on how it was built using Explore.

Andt Maple TA now supports facilities for an instructor/app-author to set certain states of embedded components within a worksheet as being gradeable. See here, perhaps.

acer

It's not clear to me whether you are trying to solve the question you asked, or some related multiobjective optimization method.

I'll try to address what you wrote in the original question, but replacing `Max` with `max`.

(Was that your only issue? That you used `Max` which means nothing to Maple...)

It looks as if `penalty` has a minimum of 0, which it attains over a triangular region. Did you have some additional weighting with which to adjust the objective, in order say to force a unique minimum?

Or were you trying to solve some other question (ie. not that `penalty`, and with g1,g2,g3 as constraints, etc)?

The result of using the objective as originally posed, after replacing Max by max, doesn't look like a sensible result if you actually wished to minimize f1 and f2 together somehow, in some multiobjective sense. Both your suggested result and the one Markiyan obtained from DirectSearch look better on that query not explicitly asked (though perhaps implied by the paper citation), with both results unsurprisingly on the upper edge of the triangle. Does the paper claim that penalty should lead to a pareto minimum?

Download penalty2.mw

acer

Can you pad the interior data points, by introducing duplicates which are only slightly off-value in the independent? That way it would technically be piecewise linear, but would almost completely behave like piecewise (near-)constant.

For example,

x:=[0.0,3.0,6.0,9.0]:
y:=[4.0,5.0,7.0,11.0]:

plots:-pointplot(zip(`[]`,x,y),style=line,view=0..11,scaling=constrained);

padx:=[x[1],seq([x[i]-1e-9,x[i]][],i=2..nops(x))]:
pady:=[seq([y[i],y[i]][],i=1..nops(x)-1),y[-1]]:

plots:-pointplot(zip(`[]`,padx,pady),style=line,view=0..11,scaling=constrained);

You might also make each new entry in pady be ever so slightly different that its (x-wise farther) neighbor, if say the piecewise linear interpolator doesn't like zero slope. This might not be necessary. Ie,

pady:=[seq([y[i],y[i]+1e-9][],i=1..nops(x)-1),y[-1]];

You could do,

plots:-listplot(result);

acer

That unknown _B1~ has the assumption that it should be either 0 or 1.

Substituting either of those two assumed values of _B1~ into your test gives something which simplifies to a trivial identity in both cases.

restart:

eq := sin(x^2)=1/2:
_EnvAllSolutions := true:

s := solve(eq);
                           1     (1/2)    1     (1/2)   
                      s := - (%1)     , - - (%1)        
                           6              6             
                    %1 := 24 Pi _B1~ + 72 Pi _Z1~ + 6 Pi

test := subs(x=s[1],eq);

                           /2                       1   \   1
                test := sin|- Pi _B1~ + 2 Pi _Z1~ + - Pi| = -
                           \3                       6   /   2

vars:=[ op( indets(test,name) minus {constants} ) ];

                            vars := [_B1~, _Z1~]

about(vars[1]);

  Originally _B1, renamed _B1~:
    is assumed to be: OrProp(0,1)


expand( subs(vars[1]=0, test) );

                                    1   1
                                    - = -
                                    2   2

expand( subs(vars[1]=1, test) );

                                    1   1
                                    - = -
                                    2   2

acer

If you don't upload the worksheet we can only guess, with so few details to go on.

There might be code in the Startup Code region.

There might be code in a collapsed Code Edit Region.

There might be action code behind an embedded component.

There might be code in a collapsed Execution Group within a Document Block.

And within or without any of those there might be calls to procedures or modules stored in separate .mla Library archive files.

acer

As Axel suggested, the imaginary component is zero.

Following that idea, one can simplify the expression (RHS-LHS of the equation), and justifiably consider its real component only. There appears to be a root where the expression's purely real value passes (continuously) through zero near dd = -2.98 or so.

restart;

TTot := 70:
TC := 17:
GM := .26:
QMax := 870:
V := 3600*GM*QMax*TTot:

Eq := V = Int(QMax*exp((-t+TC)/dd)*(1+(t-TC)/TC)^(TC/dd), t = 0 .. TTot):

# The imaginary part is zero
simplify(evalc(Im((rhs-lhs)(Eq))));

                               0

eq := value(Eq):
expr:=Re( simplify((convert((rhs-lhs)(eq),rational)),size) ):

H:=proc(x)
  if not type(x,numeric) then
    return 'procname'(args);
  end if;
  evalf(eval(expr,dd=x));
end proc:

Digits:=20:

sol := fsolve(Re(H(x)), x=-10..0);

                     -2.9847453623107653100

#plot(H(x),x=-9..-2);

evalf[1000](subs(dd=sol,Re(eq))):
evalf(%);

                        7                           7
          5.700240000 10  = 5.7002399999999999996 10 

The expression does not exceed about -5.7e7 for dd>0.

acer

Your expressions have multiple (at least two) different assumed names that print as p~.

indets(Q[6][2,3],name);

                           {p~, p~, x}

Those do not combine or cancel under `simplify` (which is usual).

I have not yet looked to see where they arise. It may be that you could use `additionally` rather than repeated calls to `assume`. Or you might be able to proceed with something crude like,

simplify(convert(Q[6][2,3],`global`));

[update]

You just have to move the line which calls `assume(p,real)` out of the P_chain procedure. It can be called just once at the top-level (since your P_chain already declares p as global), at before the loop. Here is an edited version,

June_18_modif.mw

acer

With Maple 18 one could try and use the new InertForm package for this.

The dot appearing in the imaginary term looks a little heavy to me. (It shows in the Std GUI, but not the TTY).

One nice thing is that the extra brackets don't appear, when done in the GUI.

For example,

  restart:

  c:=a+b*I:
  expr:=c^2:

  U := `%+`(evalc(Re(expr)),evalc(I*Im(expr))):

  with(InertForm):

  Display(U);

  Display(U, inert=false);

  Value(U);

acer

Are you using deprecated matrix rather than Matrix, and if so why?

Are you using Maple 18 and have you declared D as a top level local? If you are using name D at the top level (ie. this code is not appearing inside any procedure of your own) then perhaps it would simplify your issue if you switched that to say DD. 

How about something simpler such as this, if the Q[i] are all Matrices,

DD := add( evalf( (1/3)*Q[i]*(z(i+1)^3-z(i)^3) ), i=1..4 );

Note that even if you kept the loop, rather than use add, then you would not need evalm if you used Matrix instead of matrix.

If you still get problems with the diagonal elements, then perhaps you code upload the code so that we could investigate directly.

acer

I don't yet see a way to get an exact chatcterization of all the roots.

You may be able to obtain float approximations to roots. It may help to scale the expression (according to its magnitude on the requested range? Unclear to me.)

restart:

expr:=-9999990000000000000000*cos(166*beta)*sinh(166*beta)*cosh(88*beta)^2
 -9999990000000000000000*cos(88*beta)^2*sin(166*beta)*cosh(166*beta)
 +9999990000000000000000*sinh(166*beta)*cos(88*beta)^2*cos(166*beta)
 +9999990000000000000000*cosh(88*beta)^2*cosh(166*beta)*sin(166*beta)
 +10000010000000000000000*cos(166*beta)*sinh(166*beta)
 +10000010000000000000000*sin(166*beta)*cosh(166*beta)
 +9999990000000000000000*sinh(88*beta)*cos(166*beta)*cosh(166*beta)*cosh(88*beta)
 -9999990000000000000000*sinh(88*beta)*sin(166*beta)*sinh(166*beta)*cosh(88*beta)
 +9999990000000000000000*sin(88*beta)*cos(88*beta)*sinh(166*beta)*sin(166*beta)
 +9999990000000000000000*cos(88*beta)*cos(166*beta)*sin(88*beta)*cosh(166*beta)
 -9980010000000000000000*cosh(88*beta)^2*sinh(166*beta)*cos(88*beta)^2*cos(166*beta)
 -9980010000000000000000*cosh(88*beta)^2*cosh(166*beta)*sin(166*beta)
 *cos(88*beta)^2
 +9980010000000000000000*sinh(88*beta)*cos(88*beta)^2*sin(166*beta)
 *sinh(166*beta)*cosh(88*beta)
 -9980010000000000000000*cos(88*beta)*cosh(88*beta)^2*sin(88*beta)
 *sin(166*beta)*sinh(166*beta)
 +9980010000000000000000*sinh(88*beta)*cosh(88*beta)*cosh(166*beta)
 *cos(88*beta)^2*cos(166*beta)
 +9980010000000000000000*cosh(88*beta)^2*cosh(166*beta)*cos(88*beta)
 *sin(88*beta)*cos(166*beta)-9980010000000000000000*cos(88*beta)
 *sinh(88*beta)*cos(166*beta)*sin(88*beta)*sinh(166*beta)*cosh(88*beta)
 +9980010000000000000000*cos(88*beta)*cosh(88*beta)*sin(88*beta)
 *sin(166*beta)*cosh(166*beta)*sinh(88*beta):

#Digits:=500:
plot(expr,beta=0.0..0.2, view=0..1e-100);
#Digits:=10:

findroots:=proc(expr,a,b,{guard::posint:=5,maxtries::posint:=50})
local F,x,sols,i,res,start,t;
   x:=indets(expr,name) minus {constants};
   if nops(x)>1 then error "too many indeterminates"; end if;
   F:=subs(__F=unapply(expr,x[1]),__G=guard,proc(t)
      Digits:=Digits+__G;
      __F(t);
   end proc);
   sols,i,start:=table([]),0,a;
   to maxtries do
      i:=i+1;
      res:=RootFinding:-NextZero(F,start,
                                 'maxdistance'=b-start);
      if type(res,numeric) then
         sols[i]:=fnormal(res);
         if sols[i]=sols[i-1] then
            start:=sols[i]+1.0*10^(-Digits);
            i:=i-1;
         else
            start:=sols[i];
         end if;
      else
         break;
      end if;
   end do;
   op({entries(sols,'nolist')});
end proc:

findroots(expr, 0, 0.2); # missing some hinted at by the plot

  0.01641624893, 0.03039949207, 0.05572565081, 0.06691672633, 0.1327837617


findroots(expr*1e-90, 0, 0.2); # agreeing with the plot

 0.01641624893, 0.03039949207, 0.05572565081, 0.06691672633, 0.09465903392, 

   0.1039052938, 0.1327837505, 0.1327837510, 0.1417059480, 0.1700140953, 

   0.1804062466

acer

expr1 := p+p^(-1/(theta-1))*sum(q[i]^(theta/(theta-1)),i=(1..n));

                                         /  n                  \
                           /      1    \ |-----     /  theta  \|
                           |- ---------| | \        |---------||
                           \  theta - 1/ |  )       \theta - 1/|
             expr1 := p + p              | /    q[i]           |
                                         |-----                |
                                         \i = 1                /

expr2 := p^(-1/(theta-1))*(p^(theta/(theta-1))
         +sum(q[i]^(theta/(theta-1)),i=1..n));

                              /               /  n                  \\
                /      1    \ | /  theta  \   |-----     /  theta  \||
                |- ---------| | |---------|   | \        |---------|||
                \  theta - 1/ | \theta - 1/   |  )       \theta - 1/||
      expr2 := p              |p            + | /    q[i]           ||
                              |               |-----                ||
                              \               \i = 1                //

H:=(ee,a) -> a*simplify(1/a*ee):

H(expr1, p^(-1/(theta-1)));

                          /               /  n                  \\
            /      1    \ | /  theta  \   |-----     /  theta  \||
            |- ---------| | |---------|   | \        |---------|||
            \  theta - 1/ | \theta - 1/   |  )       \theta - 1/||
           p              |p            + | /    q[i]           ||
                          |               |-----                ||
                          \               \i = 1                //

Did you want to also programmaticaly "discover" that second argument to H above, the special term, perhaps as the coefficient of the Sum subexpression?

acer