acer

32333 Reputation

29 Badges

19 years, 322 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

So, are you just wondering about the formats for floats that printf supports?

If all your output intended for export is going to be produced in such a way that you might use fprintf rather than printf and writeto.

The printf command accepts HFloats just as it accepts "software floats", as input for its supported float formats. No need to convert it first, if I am understanding your question.

For example,

x:=HFloat(2.25688889345634):

printf("x=%.3f\n", x);
x=2.257

printf("x=%.3e\n", x);
x=2.257e+000

printf("x=%.3g\n", x);
x=2.26

printf("x=%Y\n", x);
x=40020E1BC3A0CE0C

acer

If you wish for hue or more easily customizable piecewise-linearly-interpolated shading schemes then you might want to consider the 2-dimensional output from plots:-surfdata that are possible with Maple 18. That is the result R in the code below.

It also gives a more smoothly interpolated default color rendering (for a hue scheme, at least).

Apart from comments in the code below, I could also mention that plots:-densityplot might be also be used to get a hue shading scheme (possibly by constructing a procedure which floored x and y arguments and returned the appropriate Matrix entry values).

And, unlike plots:-listdensityplot, the plots:-densityplot command allows more control over the end-values via its restricttoranges option. But this is also easily handled by including or omitting an offset as in the code below for calling plots:-surfdata.

And you can also quite easily make plots:-surfdata use your own customized piecewise linearly interpolated color scheme in Maple 18. So alternatively you could easily have it shade from orange through white to green, etc. (You can get plots:-densityplot to produce the same effect, but it's more effort to set up.)

NB. The 3D surfdata plot looks much better in the Maple GUI than it does inlined below. But even in the GUI the hue scheme is not a good match at one extreme. And the 2D surfdata plot also looks cleaner and smoother in the Maple 18 Std GUI than it does inlined below.


restart:

with(plots): with(plottools):
#plots:-setoptions(size=[250,250]);

N:=50:

M:=Matrix(N,N,(i,j)->evalhf(sin(i*2*Pi/N)*(j-N/2)), datatype=float[8]):


# By default, listdensityplot has tickmarks aligned with the middle
# or Matrix column- or row-values

A:=listdensityplot(M, colorstyle=HUE, thickness=0, style=surface):
A;


# 2D coordinate transformation of above plot A

f:=transform( (x,y)->[2*x+1,y-3] ):
T:=display(f(A), axes=frame, gridlines=false):
T;

#display(Array([A, T]));


# 3-dimensionsional, with forced orientation.
# Note the 3D zhue shading goes from magenta to red, thus not a true match!
# Offsets of 0.5 are used, to get same mid-point alignment.
# Note the color interpolation during rendering, smoother than above.

a,b,c,d:=(0.5)*2+1,(N+0.5)*2+1,(0.5)-3,(N+0.5)-3:

CL:=plots:-surfdata(-M, a..b, c..d, shading=zhue, style=patchnogrid,
                    orientation=[-90,0,0], lightmodel=none):
CL;


# 2-dimensionsional (Maple 18), producing a densityplot structure,
# but with finer control over the color shading.
# Note that this scheme hits red at both extremes, 0 and 1, matching
# listdensityplot HUE shading.
# Same offsets of 0.5 as above are used.
# Note the smoother color interpolation during rendering.

a,b,c,d:=(0.5)*2+1,(N+0.5)*2+1,(0.5)-3,(N+0.5)-3:

col:=["red","magenta","blue","cyan","green","yellow","red"]:
m:=[seq(i,i=0..nops(col)-1)]/(nops(col)-1.0):

R:=plots:-surfdata(-M, a..b, c..d, dimension=2, axes=frame, style=patchnogrid,
                   colorscheme=["zgradient",col,markers=m]):
R;

#display(Array([T, R, CL]));

 


Download listdens.mw

acer

Below I replace sum by Sum, scale up h by a factor of 10^1560, and increase the working precision. A result near 15.7 then comes out pretty quickly.

restart:

f := (a, b, k) -> (Sum(exp(-x)*x^m/(m!), m = a .. b))^k:

h := f(0,0,82)*f(1,3,49)*f(4,6,47)*f(7,10,47)*f(11,15,57)*f(16,20,40)
     *f(21,25,38)*f(26,35,52)*(1-f(0,35,1))^91:

H:=proc(X,d)
     if not type(X,numeric) then
       return 'procname'(args);
     end if;
     Digits:=d;
     evalf(subs(x=X,10^(1560)*h)):
end proc:

#plot(H(x,2000), x=14..18);

Optimization:-Maximize(H(x,2000),x=14..18);

               [                   5                        ]
               [6.96925226027268 10 , [x = 15.6881165583629]]

acer

The problem is your use of `sum`.

The `plot` command is receiving the result of subs(M=11,f(t)) which is expensive to compute (it is not evalhf'able). Issue the command,

subs(M=11,f(t))

to see what you were actually passing to `plot`.

So you have been hit by premature evaluation, where f(t) is evaluated as normal as part of an expression passed as the first argument to `plot`.

The problem for this example is not that Maple attempts symbolic summation for each new numeric value of `t` to be plotted (which Axel may have surmised). That situation is also quite common, and the symptoms and fix are often the same. But, technically, in this example the problem is that the initial symbolic summation succeeds and produces an expression which subsequently evaluates expensively at numeric values of `t`.

There are several alternatives, including,

# Use the inert `Sum`
restart:                                         
f:=t->Sum(-4/(Pi*(2*n+1))*sin((2*n+1)*t),n=0..M):
plot(subs(M=11,f(t)),t=0..2*Pi);

# Pass M as an extra argument to `f`.
# Here single (uneval) right-quotes around the call to `f`
# delay evaluation.
restart:
f:=(t,M)->add(-4/(Pi*(2*n+1))*sin((2*n+1)*t),n=0..M):
plot('f'(t,11),t=0..2*Pi);

# A variant of the above (fast as long as numeric M is "small enough" for
# `sum` to dispatch to `add` rather than try symbolic summation.
# Here again quotes delay the evaluation of `f`.
restart:                                             
f:=(t,M)->sum(-4/(Pi*(2*n+1))*sin((2*n+1)*t),n=0..M):
plot('f'(t,11),t=0..2*Pi);

# A variant of the first above, with `f` accepting 2 args and not using `subs`.
# Since the inert form `Sum` is used evaluation needn't be delayed.
# This is the one that Axel showed.
restart:
f:=(t,M)->Sum(-4/(Pi*(2*n+1))*sin((2*n+1)*t),n=0..M):
plot(f(t,11),t=0..2*Pi); 

acer

Do you want B_interp(N,M) to returned unevaluated if N and M are not both numeric?

B_interp := proc(W,T)
   if not ( type(W,numeric) and type(T,numeric) ) then
      return 'procname'(args);
   else
      CurveFitting:-ArrayInterpolation([FC_Map_W,FC_Map_T],FC_Map,
                                       Array(1 .. 1, 1 .. 1, 1 .. 2, [[[W, T]]]),
                                       method=linear);
   end if;
end proc:

In such a way, you should be able to construction an expression containing B_interp(N,M) where N and M are still symbols. But you could later evaluate that expression with numeric values for M and N. Eg,

  expr := ... + B_interp(N,M) + ...
  eval( expr, [N=5.4, M=3.2,...] );

Or you could `unapply` your expression with respect to M,N,etc.

acer

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?

 

restart:

with(Optimization):
f1 := -2*x1 - x2:

f2 := -x1 - 4*x2:

g1 := 2*x1 + 3*x2 - 6:

g2 := -x1:

g3 := -x2:

penalty := lambda1*max(f1-M,0) + lambda2*max(f2-M,0)
           + (M^2)*(max(g1,0) + max(g2,0) + max(g3,0)):

k := 1:

s := 1:

obj := eval(penalty,[lambda2=0.5,lambda1=0.5,M=1]);

ans := Minimize(obj,x1=-10..9,x2=-10..9);

       obj := 0.5 max(0, -2 x1 - x2 - 1) + 0.5 max(0, -x1 - 4 x2 - 1)

          + max(0, 2 x1 + 3 x2 - 6) + max(0, -x1) + max(0, -x2)

        ans := [0., [x1 = 0.266505412954423, x2 = 0.573107578136192]]

ans[2];

              [x1 = 0.266505412954423, x2 = 0.573107578136192]

eval(f1,ans[2]), eval(f2,ans[2]);

                    -1.10611840404504, -2.55893572549919

## There are several ways to find see region where `penalty`=0 (its minimum).

#plot3d(obj,x1=-1..4,x2=-1..2);

#plots:-implicitplot(obj,x1=-1..4,x2=-1..3,
#                    gridrefine=3,view=[-1..4,-1..3]);

sol := solve(convert(obj,rational));

                   /                                   2       \
           sol := { 0 <= x1, 0 <= x2, x1 <= 3, x2 <= - - x1 + 2 }
                   \                                   3       /

# Your suggested value also lies on the edge of that region.
eval( x2<=-2/3*x1+2, [x1 = 1.551123, x2 = 0.965918] );

                           0.965918 <= 0.965918000

plots:-display(
  plots:-inequal(sol,x1=-1..4,x2=-1..3,color=blue),
  plots:-pointplot([[1.551123, 0.965918]],symbol=solidcircle,
                     symbolsize=15,color=red),
  gridlines=false);

 

 

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

First 238 239 240 241 242 243 244 Last Page 240 of 336