pagan

5147 Reputation

23 Badges

17 years, 23 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are answers submitted by pagan

You might also be interested in the syntax for producing an instance of the original Array in which some particular column or row has been replaced.

Using Markiyan's data,

B := Array(1..2,1..3,[[1, 2, 3],[4., 5., 6.4]]);

                                  [ 1   2    3]
                             B := [           ]
                                  [4.  5.  6.4]

f:=x->x^2:

f~(B[2,..]);
 
                             [16., 25., 40.96]

B; # notice: unchanged

                                [ 1   2    3]
                                [           ]
                                [4.  5.  6.4]

BB:=copy(B): # or, use B itself if you want it changed in-place
BB[2,..]:=f~(BB[2,..]):
BB;

                              [  1    2      3]
                              [               ]
                              [16.  25.  40.96]

As far as I can tell, using ArrayTools:-Alias does not give any easier way in which to map along a row of a Fortran_order 2D Array or Matrix (or, similarly, along a column of a C_order 2D Array or Matrix)

It might be nice if the RowOperation and ColumnOperation commands (or some new commands) could provide the functionality to easily update rows or columns of 2D rtables. For example it might be adjusted so that the current parameter `s`, now taken only as an expression, could be accepted as a procedure to apply to the relevent entries. Such rtable "iterators" could also be allowed for updating bands or diagonals.

Since in the end you're just evaluating numerically, using r=0.2, then there is not so much of a need to produce F(r) for a nonnumeric, symbolic `r` and then only afterwards evaluating at r=0.2. Also, your call to lowercase `int` with an integration range of x=0..`rob` for floating-point `rob` is going to eventually involve a numerical integration.

With that in mind, perhaps one easy way is to try doing it purely numerically. D(IN) will return unevaluated here, and evalf(D(IN)(value)) will call fdiff the numerical differentiator.

> ING := D(IN):

> evalf(ING(.2));

                                                  -14  
                    -0.6067953678 + 3.272366947 10    I

> # that small imaginary component is just due to numeric roundoff.
> # evaluating at higher working precision makes it more negligible.

> evalf[20](ING(.2));

                                                              -25  
            -0.60679536790134102510 - 5.3295064630932992395 10    I

> # we can remove the imaginary term (safely, ie. only if "small" enough)

> simplify(fnormal(evalf(ING(0.2))));

                                  -0.6067953678

If you don't want to see tha small imaginary artefact, at all, then you could try changing your definition of IN.

> IN:=r->abs(F(r))^2/N;

                                 2
                           |F(r)| 
                      r -> -------
                              N   

> ING := D(IN):

> evalf(ING(.2));
                         -0.6067953678

Another reason for using that modified version of IN is that it only computes F(r) once, but your version computed it twice (not an actual problem here, as evalf/Int may store the first result in its remember table, but maybe good to realize and be aware of for some future example).

There are other ways to get the result, but you'd have to be careful if you use `diff`, including avoiding the error of having it try to differentiate w.r.t a number.

Part of the confusion here is the way in which you form `h`, `H`, `F`, and `IN` all as operators which call the previous one. And preassigning all those parameters with floating-point values, at the top of the sheet, does not help with that approach. Those two aspects of writing worksheets are extremely popular amongst new users, but in my opinion they both make it harder to use Maple effectively.

If you really wanted to compute IN(r) for symbolic `r` then you might have an easir time with it by using expressions instead of those unnecessary operators. (That's just my opinion.) On the other hand, if you only ever want to evaluate `IN` or its derivative at some actual floating-point numbers then using operators is OK, but in that case don't do things like eval(IN(r),r=0.2) which first calls IN at unassigned symbolic `r`.

Replace all the square brackets in equn1 with round brackets. Square brackets are for list construction. Round brackets are expression delimiters.

Then

 'solve(equn1,beta[2])'

You are missing a multiplication sign between the bracketed terms (x+2) and (x-2)

You can still force the axis labels. Using your defn of x,

plot(x, t = -10*Unit(s) ..10*Unit(s),view=-10..50,labels=[t*Unit(s),'x'*Unit(m)]);

It opens ok using maple 15's Standard GUI (maplew.exe).

It looks like it might be a .mw worksheet which has been saved or uploaded with the wrong (.mws) filename extension.

So, your Windows thinks it should open it with the Classic GUI, but the Classic GUI doesn't understand it (because it is not actually.mws).

So you could just rename it, or open with the standard gui and save as .mw file.

fsolve may not have been satisfied that at the working precision it was not able to attain some tolerance. It's a flaw, that precision and tolerance are not separate for fsolve.

Using operators instead of expressions seems to use a weaker tolerance

ans:=fsolve({(lu,lv)->FDerivative2YYMain(lu,lv, 12.5)/(lv+12.5),
             (lu,lv)->FDerivative3YYYMain(lu,lv, 12.5)/(lv+12.5)},
            [0.05435690547,0.54502809258264223447]);

  [0.054356905499263310199322866010579288988739481809993666893229496592579505151375,
   0.54502809258264223446637739872970465118604413199322232443888293557355194978678]

evalf(FDerivative2YYMain(ans[],12.5));

                                   -62
                          -7.287 10  
 
evalf(FDerivative3YYYMain(ans[],12.5));

                                   -60
                         5.34464 10   

That error message often means that the code has referenced a Matrix with indices which are names which do not evaluate to integer values. For example

restart:

u:=Matrix([[3]]):

u[i,i];
Error, bad index into Matrix

The problem is that `i` was not an integer value (within the valid bounds range of just 1 to 1).

This can happen in a few ways. You might have a line of code where you intended `i` to have a value, but it actually didn't.

One common scenario is that the code calls `sum` and the first argument to the `sum` call is a formula containing things like u[i,j] say. Since `sum` does not have special evaluation rules Maple will try and evaluate u[i,j] before `i` and/or `j` the indices of summation get assigned any numeric values. And then bang. If that is your situation then try replacing calls to `sum`, for summing finitely many terms, by calls to `add` instead.

Now here is the best bit: in the Maple Standard GUI that error message is a hyperlink. If you follow its URL in a web browser then you will get to an Online help page which gives pretty much the same explanation.

The Compiler demands that local variables cannot have Array types, but your procedure references the entries of the locals E, V, and s as if they were 1-dimensional Arrays. Hence the procedure is doing something that the Compiler does not support.

Another way to say it is that local variables cannot be Arrays or Vector or Matrices, for the Compiler to accept the procedure. (And lists and sets and tables are not allowed, since Arrays, Vectors, and Matrices of hardware datatype are the indexable structures that the Compiler admits.)

If you read the help page Compiler you can see early on where it says, "Procedures that return arrays or allocate new arrays cannot be translated." That means that local variables cannot get new Array/Vector/Matrix assigned to them, inside the procedure.

The correct way to do what you are trying is to create all the Arrays outside of the procedure, and then to pass them in as arguments. This means locals V, E, s, and also your current global u.

It's then your choice whether to have the procedure continue to return s[Nr-1], or to just return NULL (since outside of procedure X you have access to the Array you pass in for argument s, and it gets updated in-place).

You should consider making the other globals into arguments too. Using globals when not necessary just leads to sloppy code with problems, IMNSHO. But also, for the Compiler it means callbacks to main Maple and conversions from software to hardware scalar types, and if that occurs inside a loop (not your situation here...) that can be costly (relative to compiled operations).

Also, if you are going to use the Compiler than give your procedure's arguments and locals the type information (using :: double-colon) the makes things unambiguous. Don't rely on the Compiler to magically guess that you intend on using this as a float and that as an integer, etc, because if deduced wrongly the code may compute wrongly through no fault of its own.

You'll have to choose appropriate ranges, eg, alpha=0.0 .. 1.0, and n=1.., etc, in the Explore pop-up window

restart:

K := proc (alpha, theta1, theta2, sigma, n::posint)
local t,beta,W;
uses Statistics;
t[alpha] := fsolve(CDF(RandomVariable(NormalDistribution(0, 1)), -t) = alpha, t);
beta := evalf(CDF(RandomVariable(NormalDistribution(0, 1)),
t[alpha]-sqrt(n)*abs(theta2-theta1)/sigma)); 
W := evalf(CDF(RandomVariable(NormalDistribution(0, 1)),
 -t[alpha]+sqrt(n)*abs(theta2-theta1)/sigma));
sprintf("alpha is %a, beta is %a, and W is %a", alpha,beta,W);
end proc:

K(0.5e-1, 0, .6, 1, 5);

  "alpha is .5e-1, beta is .6191361682, and W is .3808638318"

eval(subs(z=eval(K),'Explore(z(alpha,theta1,theta2,sigma,n))'));

The methodology within K might be improved slightly. But that is incidental.

It looks like a regression bug.

Maple 12.02, 13.02, and 14.01 (32bit and 64bit) for Windows produces all three of these plots below correctly. But in Maple 15.01 (32bit and 64bit), the first two are produced incorrectly with either all x- or all y-values as zero. (The plot data structures contain the wrong zero-values. It's not just a GUI rendering thing.)

restart:

plot(Matrix([[1e-10,1e-9],[4e-10,4e-9],[5e-10,5e-9],[6e-10,6e-9],[9e-10,9e-9]]),
     style=point,symbolsize=20,view=[0..1e-9,0..9e-9]);

Digits:=11:

plot(Matrix([[1e-10,1e-9],[4e-10,4e-9],[5e-10,5e-9],[6e-10,6e-9],[9e-10,9e-9]]),
     style=point,symbolsize=20,view=[0..1e-9,0..9e-9]);

Digits:=12:

plot(Matrix([[1e-10,1e-9],[4e-10,4e-9],[5e-10,5e-9],[6e-10,6e-9],[9e-10,9e-9]]),
     style=point,symbolsize=20,view=[0..1e-9,0..9e-9]);

The notation of c'x for the objective is just shorthand for a Vector dot product. It's a scalar, the result of the dot product of the row Vector c' (ie. the transpose of column Vector c) and the column Vector x. The entries of Vector c are the numeric coefficients, and the Vector x are the as-yet-unknown variables. Hence c'x is a linear combination -- not surprising as we are talking about Linear Programming problems.

So, in Matrix form, the Vector x for you example will have as it entries all the w[i] as well as dr and ur. But since for your example the objective c'x is just the linear combination ur + dr = 0*w[1]+0*w[2]+...+0*w[50]+1.0*ur+1.0*dr then it is easy to see how to construct both c and x as Vectors.

The Vector c might be taken as Vector(52,i->`if`(i=51 or i=52, 1.0, 0.0))

The Vector x would then be taken as Vector(52,i->`if`(i=51,dr,`if`(i=52,ur,w[i])))

c:=Vector(52,i->`if`(i=51 or i=52, 1.0, 0.0)):
x:=Vector(52,i->`if`(i=51,dr,`if`(i=52,ur,w[i]))):

c^%T . x;
                        1.0 dr + 1.0 ur

It is trivial to adjust, for a maximation problem.

And now that you know what are the entries of the Vector x of unknowns, you can see the LPSolve Matrix form help page for how to set up the linear constraints and simple bounds.

If you are interested in efficiency, then use the syntax z[i,1] to access the entries of Matrix z, as the syntax z[i][1] is quite an inefficient way to get the same thing. When you invoke z[i][1] you force Maple to waste resources extracting z[i] (which is the i'th row of z) as a new Vector object whose transient existence serves only to have its first entry extracted. But the syntax z[i,1] extracts the first entry in the i'th row directly, as a mere scalar. So, one way forms a new Vector (which takes resources to form as well as to collect and dispose of as garbage) just to extract a single scalar entry, while the other method simply extracts that same desired scalar entry directly with no such addtional overhead.

Take notice that you cannot use `I` in that manner, unless you first switch the imaginaryunit. So I use II instead.

eq:=(Ls+R+1/Cs)*II(s)=V(s);

                   /         1 \             
                   |Ls + R + --| II(s) = V(s)
                   \         Cs/             

First, one of the simpler ways,

isolate(eq,II(s));

                                 V(s)    
                      II(s) = -----------
                                       1 
                              Ls + R + --
                                       Cs

%/V(s);

                      II(s)        1     
                      ----- = -----------
                      V(s)             1 
                              Ls + R + --
                                       Cs

which together are

isolate(eq,II(s))/V(s);

                      II(s)        1     
                      ----- = -----------
                      V(s)             1 
                              Ls + R + --
                                       Cs

Also possible here is

eq/II(s);

                               1    V(s) 
                      Ls + R + -- = -----
                               Cs   II(s)

1/rhs(%)=1/lhs(%);

                      II(s)        1     
                      ----- = -----------
                      V(s)             1 
                              Ls + R + --
                                       Cs

And, if you like typing

isolate(II(s)/V(s)=T,II(s));

                         II(s) = V(s) T

subs(%,eq);

                  /         1 \              
                  |Ls + R + --| V(s) T = V(s)
                  \         Cs/              

isolate(%,T);

                                 1     
                        T = -----------
                                     1 
                            Ls + R + --
                                     Cs

subs(T=II(s)/V(s),%);

                      II(s)        1     
                      ----- = -----------
                      V(s)             1 
                              Ls + R + --
                                       Cs

which together are

subs(T=II(s)/V(s),isolate(subs(isolate(II(s)/V(s)=T,II(s)),eq),T));

                      II(s)        1     
                      ----- = -----------
                      V(s)             1 
                              Ls + R + --
                                       Cs
f:=l^2/(m^2*M*G*(1+((m^3*G^2*M^2+2*En*l^2)/(m^3*G^2*M^2))^(1/2)*cos(phi)));

                               2                       
                              l                        
        -----------------------------------------------
               /                        (1/2)         \
               |    / 3  2  2         2\              |
         2     |    |m  G  M  + 2 En l |              |
        m  M G |1 + |------------------|      cos(phi)|
               |    |      3  2  2     |              |
               \    \     m  G  M      /              /

g:=-l^2*(M*G*m^2-(m*(m^3*G^2*M^2+2*En*l^2)/cos(phi)^2)^(1/2)*cos(phi)^2)
   /(m*(2*En*l^2*cos(phi)^2-m^3*G^2*M^2+m^3*G^2*M^2*cos(phi)^2));

         /                                 (1/2)          \
         |         /  / 3  2  2         2\\               |
       2 |     2   |m \m  G  M  + 2 En l /|              2|
      l  |M G m  - |----------------------|      cos(phi) |
         |         |              2       |               |
         \         \      cos(phi)        /               /
    - -----------------------------------------------------
        /      2         2    3  2  2    3  2  2         2\
      m \2 En l  cos(phi)  - m  G  M  + m  G  M  cos(phi) /

simplify(f-g) assuming phi>0, Pi/2>phi, m>0, G>0, M>0;

                               0

# trying to get f from g
(Z->numer(Z)/(op(1,denom(Z))*expand(denom(Z)/op(1,denom(Z)))))
(1/simplify(rationalize(denom(g)/numer(g)))) assuming phi>0, Pi/2>phi, m>0, G>0, M>0;

                               2                       
                              l                        
        -----------------------------------------------
               /                                 (1/2)\
               |             / 3  2  2         2\     |
             2 |    cos(phi) \m  G  M  + 2 En l /     |
        M G m  |1 + ----------------------------------|
               |                     (3/2)            |
               \                M G m                 /

# trying to get g from f
(Z->-collect(combine(expand(-m*numer(Z))),l)/(m*denom(Z)))
(rationalize(f)) assuming m>0, G>0, M>0;

       /                               (1/2)         \     
       |           / 3  2  2         2\              |     
       |     2     |m  G  M  + 2 En l |              |  2  
       |M G m  - m |------------------|      cos(phi)| l   
       \           \         m        /              /     
    - -----------------------------------------------------
        /      2         2    3  2  2    3  2  2         2\
      m \2 En l  cos(phi)  - m  G  M  + m  G  M  cos(phi) /
First 7 8 9 10 11 12 13 Last Page 9 of 48