acer

32622 Reputation

29 Badges

20 years, 44 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

You should be able to handle this in one of two ways. You can either test that the result is of a type (that a solution would take), or you can test whether the result is NULL or an unevaluated `fsolve` call.

result := fsolve(...);

if result = NULL or type(result, 'specfunc(anything,fsolve)' then
  no_solution_found := true;
end if;

# or, in the multivariate case,

if not type(result, set(symbol=complex(numeric))) then
  no_solution_found := true;
end if;

Notice that if you are choosing to test whether the result has a given form then the type-check may depend on how you called fsolve. These next two should illustrate that. Note also that you can choose to test whether the numeric result is purely real (ie. of type numeric which would be purely real instead of type complex(numeric) which might have nonzero imaginary part).

> result := fsolve( sin(5*x+3) , {x});
                      {x = 0.02831853072}

> type(result,set(symbol=numeric));
                              true

> result := fsolve( sin(5*x+3) , x);
                         0.02831853072

> type(result,numeric);
                              true

One can see, that it's a bit easier to check whether the `fsolve` attempt failed. If you had success, then how you pick apart the result will depend on how you called fsolve, and perhaps on how many variables were involved or whether the input was a univariate polynomial (in which case multiple solutions may be returned). How you test for that, or how you pick apart the solution, will thus depend on the problem at hand.

Now here is a related digression. You have realized that the result maybe a so-call unevaluated fsolve call, which is a fancy way of describing the situation that fsolve returns something that looks like the original call -- ie. a function call to `fsolve` itself.

If you test this at the top-level, outside of a procedure, then Maple will try once again to evaluate that function call when you inspect it by passing it on to your test. You likely don't want that, since you don't want your code to waste time duplicating a non-fruitful computation.

But if you instead perform the original fsolve call and the subsequent check inside a procedure, where result is a declared local variable, then any unevaluated fsolve function call will not get re-evaluated when your pass it around for the test.

Let's verify that claim, by simply counting how many times the argument to fsolve itself gets called. We'll use a procedure `f` as the thing to be solved-for, but the overall effect is the same even if the thing to be solved-for is an expression.

> restart:

> f:=proc(x)
>   global count;
>   count := count+1;
>   sin(5*x+3)-4;
> end proc:

# First let's do it inside a procedure.
 
> p:=proc(F) > global count; > local result; > result := fsolve(F); > printf("total number of f calls done: %a", count); > if type(result,numeric) then > print("root found"); > else > print("no root found"); > end if; > printf("total number of f calls done: %a", count); > end proc: > count := 0: > p(f); total number of f calls done: 828 "no root found" total number of f calls done: 828

Notice how the total number of calls to `f` did not increase due to calling `type` against the result

# Now, let's do it at the top-level, outiside if any procedure.

> count := 0: > result := fsolve(f); result := fsolve(f) > printf("total number of f calls done: %a", count); total number of f calls done: 828 > if type(result,numeric) then > print("root found"); > else > print("no root found"); > end if; "no root found" > printf("total number of f calls done: %a", count); total number of f calls done: 1656

Look what's happened. The total count of calls to `f` has doubled, following the `type` check!

The very act of passing the unevaluated fsolve call (assigned to `result`) to `type`() at the top-level causes it to get re-evaluated and thereby duplicates the entire unproductive attempt. But inside the procedure `p`, the variable `result` is a local variable and as such gets only one-level evaluation when passed into `type`(), and so does not cause a re-computation attempt.

When done at the top-level, in the case that the computation failed and returned unevaluated, the mere act of checking it causes the failing computation to be needlessly redone. This can be quite an important distinction when one is concerned with efficiency or when the computation takes a very long time. The very same issue would arise when calling `int` to do symbolic integrals and then passing around an unevaluated `int` call, except that `int` may be safe from this effect due to its having remember tables to store prior results. But the issue can arise elsewhere, since not all Library routines which can return unevaluated have remember tables.

This is one important reason why programming in Maple (ie. writing and using procedures) is important.

acer

Or you could try the top-level `copy` command.

acer

It looks like MTM:-subs doesn't really support multiple substitutions. This seems to be the case for both sequential substitutions (more that two arguments with all last denoting substitutions), or concurrent substitutions (first argument being a set).

If you want your `with` calls to be at the top then you can use the global form :-subs() in the code portion.

Or get rid of the with(MTM) if you don't need it, as was mentioned by Kamel.

In essence, MTM:-subs is designed to behave like the Matlab command subs and not like the Maple command `subs`. They were written for different products, and they behave differently. (It's not relevent that the Matlab `subs` command is part of it's old symbolic toolbox which used Maple once upon a time. The command was still designed for Matlab use and the Mathworks made it do what they wanted and had no cause to mimic Maple's `subs` command.)

acer

I'm not sure (yet) whether this is the full solution. I didn't have the patience to see whether solve({fff},AllSolutions) would ever return a result, as it took so long.

> solve(simplify(eval({fff},{x[22]=1,x[23]=0,x[31]=0,x[32]=1,y[21]=0})),AllSolutions);

                           /     1 \ 
                          { T = --- }
                           \    _Z1/ 

> about(_Z1);

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

I got this by substituting new simple names for each of the sin & cos calls, solving, resubstituting, and trying to test whether the resubstitutions introduced inconsistency.

Oh, hang on! There is at least this larger set,

> SOL:=solve(eval({fff},[x[31]=0]),AllSolutions);

 /     1                                             \   
{ T = ---, x[22] = 1, x[23] = 0, x[32] = 1, y[21] = 0 }, 
 \    _Z1                                            /   

   /     1                                             \    /
  { T = ---, x[22] = 1, x[23] = 0, x[32] = 1, y[21] = 0 }, { 
   \    _Z2                                            /    \

          2 Pi                                                  
  T = -------------, x[22] = 1, x[23] = 0, x[32] = -1, y[21] = 0
      Pi + 2 Pi _Z3                                             

  \    /        2 Pi                                         
   }, { T = -------------, x[22] = 1, x[23] = 0, x[32] = -1, 
  /    \    Pi + 2 Pi _Z4                                    

           \    /         2 Pi                               
  y[21] = 0 },  |T = ---------------, x[22] = 0, x[23] = -1, 
           /   <     1                                       
                |    - Pi + 2 Pi _Z5                         
                \    2                                       

                      \    /         2 Pi                   
  x[32] = 1, y[21] = 0| ,  |T = ---------------, x[22] = 0, 
                       >  <     1                           
                      |    |    - Pi + 2 Pi _Z6             
                      /    \    2                           

                                  \    /          2 Pi         
  x[23] = -1, x[32] = 1, y[21] = 0| ,  |T = -----------------, 
                                   >  <       1                
                                  |    |    - - Pi + 2 Pi _Z7  
                                  /    \      2                

                                             \    /
  x[22] = 0, x[23] = 1, x[32] = -1, y[21] = 0| ,  |
                                              >  < 
                                             |    |
                                             /    \

            2 Pi                                           
  T = -----------------, x[22] = 0, x[23] = 1, x[32] = -1, 
        1                                                  
      - - Pi + 2 Pi _Z8                                    
        2                                                  

           \ 
  y[21] = 0| 
            >
           | 
           / 

> nops({SOL});
                               8

seq(simplify(eval({fff},teqs union {x[31]=0})),teqs in {SOL});
             {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}

And of course by the symmetry of Matrices `c` and `h` in `eq` we also have to check the following. (There are some solutions unique to each set, where x[31]<>0 or x[23]<>0.)

SOL2:=solve(eval({fff},[x[23]=0]),AllSolutions);
 /     1                                             \   
{ T = ---, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, 
 \    _Z9                                            /   

   /     1                                              \   
  { T = ----, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, 
   \    _Z10                                            /   

   /     1                                              \   
  { T = ----, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, 
   \    _Z11                                            /   

   /     1                                              \   
  { T = ----, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, 
   \    _Z12                                            /   

   /     1                                              \   
  { T = ----, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, 
   \    _Z13                                            /   

   /     1                                              \    /
  { T = ----, x[22] = 1, x[31] = 0, x[32] = 1, y[21] = 0 }, { 
   \    _Z14                                            /    \

           2 Pi                                                  
  T = --------------, x[22] = 1, x[31] = 0, x[32] = -1, y[21] = 0
      Pi + 2 Pi _Z15                                             

  \    /         2 Pi                                         
   }, { T = --------------, x[22] = 1, x[31] = 0, x[32] = -1, 
  /    \    Pi + 2 Pi _Z16                                    

           \    /         2 Pi                             
  y[21] = 0 }, { T = --------------, x[22] = 1, x[31] = 0, 
           /    \    Pi + 2 Pi _Z17                        

                       \    /         2 Pi                  
  x[32] = -1, y[21] = 0 }, { T = --------------, x[22] = 1, 
                       /    \    Pi + 2 Pi _Z18             

                                  \    /         2 Pi       
  x[31] = 0, x[32] = -1, y[21] = 0 }, { T = --------------, 
                                  /    \    Pi + 2 Pi _Z19  

                                             \    /
  x[22] = 1, x[31] = 0, x[32] = -1, y[21] = 0 }, { 
                                             /    \

           2 Pi                                                  
  T = --------------, x[22] = 1, x[31] = 0, x[32] = -1, y[21] = 0
      Pi + 2 Pi _Z20                                             

  \    /     1                                              \   
   }, { T = ----, x[22] = 1, x[31] = 1, x[32] = 0, y[21] = 1 }, 
  /    \    _Z21                                            /   

   /     1                                                \    /
  { T = ----, x[22] = 1, x[31] = -1, x[32] = 0, y[21] = -1 }, { 
   \    _Z22                                              /    \

           2 Pi                                                  
  T = --------------, x[22] = 1, x[31] = 1, x[32] = 0, y[21] = -1
      Pi + 2 Pi _Z23                                             

  \    /         2 Pi                                         
   }, { T = --------------, x[22] = 1, x[31] = -1, x[32] = 0, 
  /    \    Pi + 2 Pi _Z24                                    

           \    /         2 Pi                  
  y[21] = 1 }, { T = --------------, x[22] = 1, 
           /    \    Pi + 2 Pi _Z25             

                /  2    \                           /  2    \\   
  x[31] = RootOf\_Z  + 3/, x[32] = 2, y[21] = RootOf\_Z  + 3/ }, 
                                                             /   

   /     1                             /  2    \              
  { T = ----, x[22] = 1, x[31] = RootOf\_Z  + 3/, x[32] = -2, 
   \    _Z26                                                  

                 /  2    \\ 
  y[21] = -RootOf\_Z  + 3/ }
                          / 

> nops({SOL2});
                               18

> seq(simplify(eval({fff},teqs union {x[23]=0})),teqs in {SOL});
{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, 

  {0}, {0}, {0}, {0}, {0}

acer

You can either not use Units:-Standard:-`^` and keep the float powers, or you can use rational powers.

For the first of those, note that you can get all of Units:-Standard loaded excepting some individual parts such as its `^` export. For example,

> restart:
> with(Units[Standard]): unwith(Units:-Standard,`^`):

> n:=1.38:
> p1:=1.1*Unit(bar):
> v1:=3*Unit(dm^3):
> v2:=0.4*Unit(dm^3):

> v1^n;
               4.554359738 ['dm'^3]^1.38

> p2:=(p1*v1^n)/(v2^n);
               17.74096457 ['bar']

Notice how the exponents in v1^n do not simplify above. I doubt that one could get them to combine using `simplify`, either.

For the second of those two approaches mentioned, you can either just enter everything youself as pure rationals (eg. 1.38 = 138/100, etc) or you can convert them with Maple (either once up front, or on-the-fly).

> restart:
> with(Units[Standard]):

> n:=1.38:
> p1:=1.1*Unit(bar):
> v1:=3*Unit(dm^3):
> v2:=0.4*Unit(dm^3):

> R:=t->convert(t,rational):

> evalf( (R(p1)*v1^R(n))/(R(v2)^R(n)) );
                             6                  
               1.774096458 10  ['Pa']

acer

You don't have MYPLOT as the return value of your procedure. So the plotobj argument doesn't get displayed or "printed".

You can work around that in a few ways. One way might be to assign MYPLOT to a local, then reset the 'default' plot options, and then to return that local as the last line of the procedure.

Or you could do it like this,

resizeMYPLOT := proc(plotObj)
   plotsetup('jpeg', 'plotoutput' = "c://temp/bar.jpg",
             'plotoptions' = "height=1000,width=1000");
   print(plotObj);
   plotsetup('default');
   NULL;
end proc:

plot(cos(x),x=-12..12);
resizeMYPLOT(%);

acer

I tried copying from and pasting to various Components' properties fields successfully on 64bit Maple for Windows 7 using Ctl-x, Ctl-c, and Ctl-v. Those keyboard shortcuts worked for me.

acer

I wasn't sure which it is of mol or m^3 that you want to get automatically changed into NM3. Below, I'll take it as mol which you want handled that way. Presumably you want anything which would otherwise return involving mol to instead return involving NM3. Is that right?

First, one can add a new unit. After that, then you can do things like convert(..,units,mol/sec,NM3/sec).

But you can also add a new system, which is a version of SI which uses NM3 instead of mol for the dimension of `amount_of_substance`.

> restart

> Units:-AddUnit(NM3,conversion=22.413*1000*mol,context=SI,prefix=SI);

> convert(20, units, mol, NM3);

                        0.0008923392674

> Units:-AddSystem(MySI,Units:-GetSystem(SI),NM3);
> Units:-UseSystem(MySI):

> with(Units[Natural]):

> 15*10^5*mol/sec;

                           [[  NM3   ]]
               66.92544506 [[ ------ ]]
                           [[ second ]]

If I may, I advise using Units[Standard] instead of Units[Natural]. I believe that the latter is a design mistake, basically robbing one of most all single letters (which otherwise come in useful for other purposes). It can also lead to confusion. Note that Units[Standard] can be used without need for having the name of the Unit() constructor appear in 2D Math input -- palette entry or command-completion can directly get the double-bracket unit appearance for 2D Math input.

acer

There is no easy way to get such fine grained display of intermediary steps for all manner of calculations.

But if you are able to do the requisite programming then you can sometimes coerce Maple into emulating such an environment.

Fo example see here or here.

acer

Quite a few requests have been made, over several years, for some quick, easy visual cue whether an indexed name is a table reference or an atomic identifier underneath.

Would it be useful if I were to try to put together a quick context-menu entry (automatically loadable into your session using user-owned initialization files, say) which could immediately show the situation upon right-click (but without  necessitating actual 1D conversion)?

acer

The ss1 object is a Maple module. And it has exports a, b, c, and d. (See the ?exports help-page.)

Hence you can access them programmatically like so,

ss1:-a;

ss1[a];

Of the two forms above, the first works the same even if you have assigned to the global name a. But if you have assigned to global a, then the second form needs to be called like ss1['a'] instead.

And if you are inside a procedure with a local a, and global a has also been assigned, then the first form works unchanged. But the second form would have to become something like ss1[':-a'] instead.

For these reasons it is often easier to just use ss1:-a, to avoid having to type out all of ss1[':-a'] which looks silly in comparison. One of the few benefits of the ss1[a] form is that you can do nifty stuff like seq(ss1[x], x in [a,b,c,d]).

acer

@brian abraham One of the common goals of using a Window Function is to correct the oscillation phenomenon near the end-points. Since there is (by defn) no original data beyond what is supplied then the FFT will usually not be able to accomodate such a discontinuity. The FFT process relies on "knowing" that the data extends in some (overlay of) a periodic way. A hard nonzero boundary for the data thus induces an effect similar to the Gibbs phenonmenon.

A key purpose of using a Window Function to scale the data so that it tends strongly to zero near the original data point boundary is to allow the data to be be safely padded further on with zero values. Since the scaled data is converging to zero, then padding with zero (or near zero) values will make the FFT process see a natural continuation. This makes the FFT and Inverse FFT behave much better at the "old" boundary, since they now act on a wider data set. And it's not just that the "old" boundaries are well inside the extended data set. The fact that the data is all nice is also key. The augmented data can easily be made to appear to extend the old curve that is tending to zero near the old boundary. The augmented data doesn't even have to oscillate in order to effectively do its job. There is no significant distinction from the augmented data and the (now tiny) oscillation in the scaled old data near the old boundary. The augmented data is effectively periodic (though literally it is not) because it so close to zero and its would-be oscillation is miniscule.

I will try to demonstrate using the previously posted example. I've added yet another, third, data set. It is based on the Windowed data set, but padded by ten percent more points. It's not padded with zero. It's padded with a continuation of the scaled curve. I've added two more plots, to show close-up where the padded data meets the scaled original data.



 

restart:

 

freqpass:=proc(M::Vector,rng::range(nonnegint))

local Mt, Xlow, Xhi;

  Xlow,Xhi := op(rng);

  Mt := DiscreteTransforms:-FourierTransform(M):

    if Xlow>1 then ArrayTools:-Fill(Xlow-1,0.0,Mt,0,1); end if;

    if Xhi<op(1,Mt) then ArrayTools:-Fill(op(1,Mt)-Xhi-1,0.0,Mt,Xhi+1,1); end if;

    Mt[1]:=Mt[1]/2;

    2*map(Re,DiscreteTransforms:-InverseFourierTransform(Mt));

end proc:

 

f := x -> 5*x + sin(128*x) - sin(130*x) + sin(133*x) - sin(137*x);

(1)

a,b,n := -Pi,Pi,1000:

 

M := Vector(n,(i)->(evalhf@f)(a+i*(b-a)/n),datatype=float[8]):

Mg := Vector(n,(i)->evalhf(f(a+i*(b-a)/n)*exp(-1/2*(a+i*(b-a)/n)^2)),datatype=float[8]):
numpad := trunc(n*0.10):
Mgpadded := Vector(n+2*numpad,datatype=float[8]):
Mgpadded[numpad+1..numpad+1+n]:=Mg:
Mgpadded[1..numpad] := Vector(numpad,(i)->M[1]
           *evalhf(exp(-1/2*(a+(-numpad+i)*(b-a)/n)^2)),datatype=float[8]):
Mgpadded[n+numpad+1..n+2*numpad] := Vector(numpad,(i)->M[n]
           *evalhf(exp(-1/2*(a+(n+i)*(b-a)/n)^2)),datatype=float[8]):

 

plots:-pointplot([seq([k,M[k]],k=1..op(1,M))],style=line,view=[0..1000,-20..20]);

plots:-pointplot([seq([k,Mg[k]],k=1..op(1,Mg))],style=line);
plots:-pointplot([seq([k,Mgpadded[k]],k=1..op(1,Mgpadded))],style=line);
plots:-pointplot([seq([k,Mgpadded[k]],k=1..numpad+100)],style=line);
plots:-pointplot([seq([k,Mgpadded[k]],k=n+numpad+1-100..n+2*numpad)],style=line);

 

 

 

 

 

filtered := freqpass(M,1..30):

filteredg := freqpass(Mg,1..30):
filteredgpadded := freqpass(Mgpadded,1..30)[numpad+1..numpad+1+n]:

plots:-pointplot([seq([k,filtered[k]],k=1..op(1,filtered))],style=line);

plots:-pointplot([seq([k,exp(1/2*(a+k*(b-a)/n)^2)*filteredg[k]],
                      k=1..op(1,filteredg))],style=line);
plots:-pointplot([seq([k,exp(1/2*(a+k*(b-a)/n)^2)*filteredgpadded[k]],
                 k=1..op(1,filteredgpadded))],style=line);

 

 

 

 

 



Download window_correction.mw

I liked Axel's most, primarily because it gives a smooth looking curve quickly.

But just for fun, here is yet another way to get a less smooth looking curve more slowly. (I like the space that shows between the curve and the y-axis. Sigh.)

ee := evalc( subs(z=x+y*I, abs(z-2)=2 ) );
plots:-implicitplot(ee, x=-4..4, y=-4..4);

acer

Can you use Maple to demonstrate that "A_10x10 is not PD under the assumption that B_4x4 is not PD" much faster than you can demonstrate that "B_4x4 is PD under the assumption that A_10x10 is PD"?

The latter of those two ways entails computing the whole logical formula for whether A is PD, up front. You've discovered that unless is is sparse or special then this can be very expensive.

But what if the logical formula for whether B_4x4 was PD was more cheaply attainable?

And what if you also devised a test of PD which could utilize a given logical formula T at any given internal logical step, by demanding that any subcheck be done under the assumption of T?

Under such circumstances, you could use the negation of the logical formula for B_4x4 being PD as T. And you could utilize ~T while computing whether A_10x10 was PD or not. If you got lucky then that computation would return 'false' quickly, reaching a logical contradiction before having to complete. That would allow you to conclude your hypothesis.

Now, how can one compute PD of a Matrix under some logical condition? If you look at (line numbers for Maple 14),

showstat(LinearAlgebra:-LA_Main:-IsDefinite,24..38);

then you can see code for testing PD of a symmetric Matrix, I think. Near the bottom there is a test, 0 < Re(de). You might be able to change that to something like (0<Re(de) and not(T)) or similar. Or maybe it would be better to use `is` and `assuming`?

Oh, but now I realize that you are interested in positive-semidefinite and not positive-definite. I doubt that the above approach is suitable for any method that has to compute whether A_10x10 is PSD by calculating all its eigenvalues at once. Now, what was that rule for computing PSD using minors? Did it involve all the minors, as opposed to just the principal minors, or am I misremembering?

acer

You are mistaken, in your claim that "In NLPSolve, f(Pi) or f(-Pi) are also calculated". If you trace through the calculuation (using `trace`, or `stopat`, or just with `print(x)` at the start of `f`) then you can see that only floating point approximations are passed to `f` during your call to NLPSolve.

acer

First 287 288 289 290 291 292 293 Last Page 289 of 339