restart:

F:=proc(x) local y;
     y:=sprintf("%a",op(1,x));
     [seq(parse(y[i]),i=1..length(y))];
end proc:

F(0.567);
                           [5, 6, 7]

x:=evalf[100000](5555/7):

CodeTools:-Usage( F(x) ):
memory used=210.69MiB, alloc change=0 bytes, cpu time=765.00ms, real time=762.00ms

For reasons I don't fully understand, above approximately 10^5 or 10^6 digits it seems that garbage collection dominates the computation and the timing grows worse than linearly w.r.t. digits. This seems to affect both Carl's original approach and mine. I'm not sure where the culprit is, but I suspect the very many calls to builtin `parse`. I got no better with sscanf.

acer

The procedure that produces the list from the set has to be told how high it should check, no?

restart:                                          

m:=LL->{seq(`if`(LL[i]=1,i,NULL),i=1..nops(LL))}: 

unm:=(SS,k)->[seq(`if`(member(i,SS),1,0),i=1..k)]:

L:=[1,0,1,1,0,0]:                                 

S:=m(L);                                          

                                S := {1, 3, 4}

unm(S,nops(L));                                   

                              [1, 0, 1, 1, 0, 0]

acer

I don't know if there is a way to truly clear the 'value' of an already inserted PlotComponent.

I usually just code something to make the component's 'value' be something effectively close to an empty plot, using the `plot` command and without having to build a PLOT structure manually.

Eg, if the identity of the PlotComponent is "Plot0",

  DocumentTools:-SetProperty("Plot0",':-value',plot([[undefined,undefined]],':-axes'=':-none'));

So try that as the Action code behind your "Clear" Button.

acer

It's not very difficult to find two real roots of the system using `fsolve` if eq[1] is first solved for `a` and then substituted into eq[2], or if some special range for variable x is supplied.

restart:
f := x*exp((1 - x)*a):
f[2] := subs(x=f,f):
eq[1] := f[2] - x:
eq[2] := diff(f[2],x) + 1:
eq[4] := a=solve(eq[1],a);
                               /  x - 2\
                             ln|- -----|
                               \    x  /
                       a = - -----------
                               -1 + x   
R3 := subs( eq[4], eq[2] ):

#plot(R3, x=0..2);

rsols := [ fsolve(R3, x=0..1), fsolve(R3, x=1..2) ];
                  [0.2777041405, 1.722295859]

eval( subs(eq[4],[eq[1],eq[2]]), x=rsols[1] );
                      [     -10        -9]
                      [-5 10   , 1.9 10  ]

[x=rsols[1], eval( eq[4], x=rsols[1])];
              [x = 0.2777041405, a = 2.526467726]

[x=rsols[2], eval( eq[4], x=rsols[2])];
               [x = 1.722295859, a = 2.526467725]

Digits:=100:

rsols := [ fsolve(R3, x=0..1), fsolve(R3, x=1..2) ]:

eval( subs(eq[4],[eq[1],eq[2]]), x=rsols[1] );
                      [    -100       -99]
                      [1 10    , -1 10   ]

It is interesting that `fsolve` has difficulty solving the system unless it gets help with at least one of the ranges.

fsolve( [eq[1],eq[2]], {a,x}, x=0..1 );
              {a = 2.526467726, x = 0.2777041405}

fsolve( [eq[1],eq[2]], {a,x}, x=1..2 );
               {a = 2.526467726, x = 1.722295859}

And what is this nonsense? (Happens in M17.02 with Digits=10 but not Digits>10 it seems.)

s:=fsolve( {eq[1],eq[2]}, {a,x}, x=-3..0 );
              {a = 6.564724757, x = -2.865536950}

eval( [eq[1],eq[2]], s );
   [               856928220820                856928220833]
   [-1.472743519 10            , 2.009082419 10            ]

acer

acer

One of the advantages of ArrayInterpolation is that is can be used to compute interpolatory results for each input value in an Array in just a single call. If you have many values at which to interpolate, and if the input values do not depend upon each other, then this can improve performance as the problem size increases.

In your example the values at which to interpolate are the first column of Matrix `A`, and they are all known in advance of the need to interpolate. So you can bring a vectorized call to ArrayInterpolation outside of your do-loop.

InterProc1 := proc(N)
uses CurveFitting;
local i, A, variable1, variable2;
     A := Matrix(N, 4);
     variable1 := evalf([0, 100]);
     variable2 := evalf([12, 20]);
     for i from 1 to N do
          A(i, 1) := evalf(3+i);
          A(i, 2) := 2*i+A(i, 1);
          A(i, 3) := A(i, 2)-1;
          A(i, 4) := ArrayInterpolation(variable1, variable2, A(i, 1));
     end do;
     A;
end proc:

ans1 := CodeTools:-Usage( InterProc1(1000) ):
memory used=17.20MiB, alloc change=24.00MiB, cpu time=220.00ms, real time=228.00ms

InterProc1b := proc(N)
uses CurveFitting;
local i, A, variable1, variable2;
     A := Matrix(N, 4);
     variable1 := evalf([0, 100]);
     variable2 := evalf([12, 20]);
     for i from 1 to N do
          A(i, 1) := evalf(3+i);
          A(i, 2) := 2*i+A(i, 1);
          A(i, 3) := A(i, 2)-1;
     end do;
     A[1..N, 4] := ArrayInterpolation(variable1, variable2, A(1..N, 1));
     A;
end proc:

ans1b := CodeTools:-Usage( InterProc1b(1000) ):
memory used=459.64KiB, alloc change=0 bytes, cpu time=0ns, real time=11.00ms

The two results are the same.

LinearAlgebra:-Norm(ans1b-ans1);
                               0.

Some additional refinements might squeeze a bit more speed out of it (which might only be apparent for much larger sizes such as N=10^4 here). And perhaps more still could be done, reusing Vector workspaces and using ArrayTools:-Alias. But the subsequent improvements seem relatively minor compared to the improvement of making just a single call to ArrayInterpolation.

InterProc2 := proc(N)
uses CurveFitting;
local V1, i, A, variable1, variable2;
     A := Matrix(N, 4, datatype=float[8]);
     variable1 := Vector([0, 100], datatype=float[8]);
     variable2 := Vector([12, 20], datatype=float[8]);
     V1:=Vector(N, i->3+i, datatype=float[8]);
     A[1..N, 1] := V1;
     A[1..N, 2] := Vector(N, i -> 2*i+A[i,1], datatype=float[8]);
     A[1..N, 3] := Vector(N, i -> A[i,2]-1, datatype=float[8]);
     A[1..N, 4] := ArrayInterpolation(variable1, variable2, V1);
     A;
end proc:

ans2 := CodeTools:-Usage( InterProc2(1000) ):
memory used=171.74KiB, alloc change=0 bytes, cpu time=0ns, real time=7.00ms

LinearAlgebra:-Norm(ans2-ans1);
                               0.

acer

Student:-Calculus1:-Roots(cosh(C+cosh(C))/cosh(C)-2, C = -3 .. 3, numeric);

                  [-2.383208606, 0.3230741020]

acer

CurveFitting/ArrayInterpolation

acer

Your call to `animate` is passing A(tranz) which evaluates before `tranz` gets any numeric value. So `plot3d` is seeing just the valueless symbol `tranz`, and cannot produce a plot. That's what the error message is about.

This is a common usage mistake. There are several ways to correct for it. One easy way is to use so-called "uneval" quotes (single right-quotes) to delay the evaluation of the call A(tranz).

restart:

with(plots):
with(plottools):

A :=tranz-> plot3d(-2/sqrt(x^2+y^2)+5/sqrt((x-1.6)^2+y^2),
                   x = -5 .. 5, y = -5 .. 5,
                   view = [-5 .. 5, -5 .. 5, -5 .. 5],
                   transparency = tranz):
B := sphere([0, 0, 0], 2*(1/10), color = magenta, style = patchnogrid):
C := sphere([1.6, 0, 0], 5*(1/10), color = green, style = patchnogrid):
E := plot3d(0, x = -5 .. 5, y = -5 .. 5, view = [-5 .. 5, -5 .. 5, -5 .. 5],
            style = wireframe, shading = zgrayscale):

animate(display, ['A'(tranz), B, C, E, scaling = constrained,
                  view = [-5 .. 5, -5 .. 5, -5 .. 5],
                  axes = normal], tranz = 0.1 .. 0.9);

acer

rsolve( {x(n+2)-5*x(n+1)+6*x(n)=4*n+1, x(0)=1, x(1)=2}, x(n) );

                         n   7         3  n
                     -4 2  + - + 2 n + - 3 
                             2         2   

You can check that,

eval(sol,n=0);

                               1

eval(sol,n=1);

                               2

simplify( eval(sol,n=k+2)-5*eval(sol,n=k+1)+6*eval(sol,n=k) );

                            1 + 4 k

acer

You could compare the output from these various calls.

restart:
with(plots):

lambda1:=(S+sqrt(S^2+4*alpha))/(2):
lambda2:=(S-sqrt(S^2+4*alpha))/(2):

plot3d({lambda1,lambda2}, alpha=-5..5, S=-10..10);

isolate(lambda1=lambda2,alpha);

plot3d({lambda1,lambda2}, alpha=max(-1/4*S^2,-5)..5, S=-10..10);

plot3d({lambda1,lambda2}, alpha=-1/4*S^2..5, S=-10..10,
       view=[-5..5,-10..10,-10..10]);

contourplot({lambda1,lambda2}, S=-10..10, alpha=-5..5,
            contours=50, grid=[100,100]);

plottools:-transform((x,y)->[y,x])(contourplot({lambda1,lambda2},
                                                alpha=max(-1/4*S^2,-5)..5,
                                                S=-10..10, contours=50));

acer

restart:

f := x -> x^3:
g := x -> x^3:

addressof(eval(f));
                      18446744074339582038

addressof(eval(g));
                      18446744074339582126

evalb( f = g );
                             false

evalb( eval(f) = eval(g) );
                              true

is( f = g );
                             false

is( eval(f) = eval(g) );
                              true

restart:

f := x -> x^3:
g := t -> t^3:

addressof(eval(f));
                      18446744074339582038

addressof(eval(g));
                      18446744074339582126

evalb( f = g );
                             false

evalb( eval(f) = eval(g) );
                              true

is( f = g );
                             false

is( eval(f) = eval(g) );
                              true

acer

Use the currentdir command to find out where the write is being attempted. Eg, just issue

currentdir();

If that shows some location where you really shouldn't be trying to write your own files (eg. "C:/Program Files/Maple 17" then use that same command to set a more prudent location that is safer from the risk of inadvertantly clobbering one of Maple's own installed files.

For example,

currentdir(kernelopts(':-homedir'));

The location of your Maple installation is not a good place to be writing files. If you are fortunate then it will have been set write-protected by the installation process.

acer

It may happen that you want to expand the `sin` call (in the sense of using only the sum identity) but do not have at hand the terms to be temporarily frozen (via the double `subs`).

One can make a further distinction, according to which behaviour is wanted when the argument to `sin` is not a sum. It might be desired that the `sin` term be left alone, or it might be desired that it be handled in the usual `expand` manner.

Just for fun, here are two ideas. (It might be that applyrule alternatively could be used for part of it. I'd be surprised if they couldn't be improved or robustified.)

> restart:   
                                         
> expandsinsum := x -> evalindets(x,specfunc(`+`,sin),
>    q->sin(op([1,1],q))*cos(op([1,2..-1],q))
>       +cos(op([1,1],q))*sin(op([1,2..-1],q))):

> expandsinsum( 2*sin(x+Pi/4) );

                                   1/2           1/2
                           sin(x) 2    + cos(x) 2

> expandsinsum( 2*sin(2*x+Pi/4) );

                                   1/2             1/2
                         sin(2 x) 2    + cos(2 x) 2

> expandsinsum( 2*sin(2*x+4*k) );

                   2 sin(2 x) cos(4 k) + 2 cos(2 x) sin(4 k)

> expandsinsum( 2*sin(4*k) ); # not a sin of a sum; does nothing

                               2 sin(4 k)

> restart:                                                     

> expandtrigsum := x -> eval(expand( x,                        
>    op(map(op@op,indets(x,'specfunc(`+`,{sin,cos,tan})'))) )):

> expandtrigsum( 2*sin(x+Pi/4) );                              

                                   1/2           1/2
                           sin(x) 2    + cos(x) 2

> expandtrigsum( 2*sin(2*x+Pi/4) );                            

                                   1/2             1/2
                         sin(2 x) 2    + cos(2 x) 2

> expandtrigsum( 2*sin(2*x+4*k) );                             

                   2 sin(2 x) cos(4 k) + 2 cos(2 x) sin(4 k)

> expandtrigsum( 2*sin(4*k) ); # not trig of a sum; behaves like expand

                                      3
                      16 sin(k) cos(k)  - 8 sin(k) cos(k)

It's possible, then, to craft your own procedures which can produce a variety of more general behaviour.

acer

You wrote, "actually, my expectation is simple, just find back the original matrix from eigenvector and eigenvalue."

Your Matrix happens to be diagonaliable, and there are three linearly independent eigenvectors associated with its three eigenvalues. The 3x3 Matrix of its eigenvectors (as columns, like the `Eigenvectors` command returns) is invertible.

The eigen-solving, or the reconstruction, can be done for this example either as an exact symbolic or as an approximate floating-point computation.

restart:

with(LinearAlgebra):

A := `<|>`(`<,>`(1, 2,3), `<,>`(2, 3, 0), `<,>`(2, 0, 0));

                                [1  2  2]
                                [       ]
                           A := [2  3  0]
                                [       ]
                                [3  0  0]

v, EigenVector1:= Eigenvectors(A):

map(radnormal,EigenVector1.DiagonalMatrix(v).EigenVector1^(-1));

                              [1  2  2]
                              [       ]
                              [2  3  0]
                              [       ]
                              [3  0  0]

v, EigenVector1:= Eigenvectors(evalf(A)):

map(simplify@fnormal,EigenVector1.DiagonalMatrix(v).EigenVector1^(-1));

               [1.000000000  2.000000000  2.000000000]
               [                                     ]
               [         2.  3.000000000           0.]
               [                                     ]
               [3.000000000           0.           0.]

Note that not all Matrices have a full set of linearly independent eigenvectors which can form such an invertible Matrix. But even if that Matrix of eigenvectors is not invertible then the vectorized definition of eigenvalues & eigenvectors (A.x=lambda.x) can still hold. Ie, the following is very similar to the above, but both sides are not right-multiplied by the inverse of `EigenVector1` Matrix.

v, EigenVector1:= Eigenvectors(evalf(A)):

map( simplify@fnormal, A.EigenVector1 - EigenVector1.DiagonalMatrix(v) );

                            [0.  0.  0.]
                            [          ]
                            [0.  0.  0.]
                            [          ]
                            [0.  0.  0.]

v, EigenVector1:= Eigenvectors(A):

map( radnormal, A.EigenVector1 - EigenVector1.DiagonalMatrix(v) );

                              [0  0  0]
                              [       ]
                              [0  0  0]
                              [       ]
                              [0  0  0]

acer