Q: "Does it really means that there are 32950 digits between the two parts of the expression?"
A: Yes.

You can get the number of digits in a various ways. For positive integer `x` you could try length(x) or something like  ceil(log[10](x*1.0)) but you should probably test those first smaller values.

If you really want to print out all the digits then you could lineprint them (without special 2D Math output typesetting) using the `lprint` command, or if you feel brave you could try and bump up the Standard GUI's term elision threshold. But be warned that both of those (lprint, less so) can drive the Standard GUI into a nearly inoperable state. You should be able to change the term elision threshold using either the main menubar's Tools>Options (and choosing the "Precision" tab in the pop-up window) or the `interface` command.

Another possibility, for printing out a posint with bery many digits, might be to use the `writeto` command to remporarily redirect output to a plainttext file.

acer

solve(identity(a*x+b*y=2*x+5*y,x),[a,b]);

                                 [[a = 2, b = 5]]

acer

The limits as you approach A3(3) from each side are +infinity and -infinity. Which result you get, as you compute in floating-point, might depend upon the working precision (Digits) due to numerical error.

restart:
alpha_p:=1: p:=2: mu0:=4*Pi*1e-7: Br:=1.12:

A1:=unapply(sin((n*p+1)*alpha_p*Pi/(2*p))/((n*p+1)*alpha_p*Pi/(2*p)),n):

A2:=unapply(sin((n*p-1)*alpha_p*Pi/(2*p))/((n*p-1)*alpha_p*Pi/(2*p)),n):

M1:=unapply((Br/mu0)*alpha_p*(A1(n)+A2(n)),n):

M2:=unapply((Br/mu0)*alpha_p*(A1(n)-A2(n)),n):

M3:=unapply(M1(n)+n*p*M2(n),n):

A3:=unapply(((n*p-1/(n*p))*M1(n)/M3(n)+1/(n*p)),n):
evalhf(A3(3));
                                         15
                     4.595173115007179 10  
A3c:=Compiler:-Compile(A3):
A3c(3);
                                         15
                     5.514207738008614 10  
for i from 10 to 20 do
  Digits:=i;
  q:=A3(3.0):
  printf("\nDigits: %ld   evalf(A3(3)): %e\n",Digits,evalf(q));
end do:

Digits: 10   evalf(A3(3)): -2.292637e+09

Digits: 11   evalf(A3(3)): -1.604846e+11

Digits: 12   evalf(A3(3)): -8.024228e+11

Digits: 13   evalf(A3(3)): -2.292637e+12

Digits: 14   evalf(A3(3)): -2.292637e+13

Digits: 15   evalf(A3(3)): 2.674743e+14

Digits: 16   evalf(A3(3)): 1.458951e+15

Digits: 17   evalf(A3(3)): 1.604846e+16

Digits: 18   evalf(A3(3)): -4.012114e+17

Digits: 19   evalf(A3(3)): 1.783162e+18

Digits: 20   evalf(A3(3)): -3.209691e+19

A3(3);
                                     (1/2)    
                    Float(infinity) 2        1
                  - ---------------------- + -
                               2             6
                             Pi               

limit(A3(x),x=3,left);
                        Float(infinity)

limit(A3(x),x=3,right);
                        -Float(infinity)

plot(A3, 2.5..3.5);

acer

What is the value for (D@D)(f)(0) ? Is it known, as an initial condition at eta=0? Or, might you want to have it be treated as another parameter?

d2f0.mw

acer

The question asked for integer solutions so, while 0 is not allowed (due to division) and -1 is not allowed to forbid 0=2, other negative values are ok?

> restart:
> e:=(1 + 1/x)*(1 + 1/y)*(1 + 1/z) = 2:               

> simplify((rhs-lhs)(expand(eval(e,[x=1,y=-z-1])/2)));
                                       0

So for x=1 then all y=-z-1 are ok, except of course {y=0,z=-1} and {z=0,y=-1}.

In a somewhat simiular way, if a=x+1 (and a<>0 and a<>1 due to restrictions on x) then for any other integer a we want integer solutions for y and z in,

> solve((rhs-lhs)(eval(expand(e),x=a-1)),y):

> y=numer(%)/collect(denom(%),z);           
                                     a (z + 1)
                               y = -------------
                                   (a - 2) z - a

But I don't know anything better to do with that than pump in a=2,3,4,...

acer

Note that even 'linear' alone would not give phi(a+x) -> phi(a) + phi(x) if `a` is unassigned. [edit: oops, not what I intended to convey. see followup below]

Also, not that map(phi,x) evaluates to phi(x) under Maple's normal mode of evaluating arguments of a procedure call. The single right-quotes used below delay this.

> restart:

> define('phi',phi(sigma)=sigma,phi(x::{`*`,`+`})='map'(phi,x));

> phi(a*y*sigma*z);                                             

                          phi(a) phi(y) sigma phi(z)

> phi(a+y+sigma+z);

                       phi(a) + phi(y) + sigma + phi(z)

I'm not sure that I understand the behaviour you want. For unassigned `a` and `z` what do you want phi(a+z) to do? And phi(a+4*sigma+z)? [edit: oops, not what I intended to ask. see followup below]

acer

normal((3*h^2+12*h)/h);

                                           3 h + 12

acer

V:=Vector[row](8,(i)->i);

                        V := [1, 2, 3, 4, 5, 6, 7, 8]

with(ArrayTools):

m1:=Reshape(V,[2,4]);

                                   [1  3  5  7]
                             m1 := [          ]
                                   [2  4  6  8]

m2:=Alias(V,[2,4]);

                                   [1  3  5  7]
                             m2 := [          ]
                                   [2  4  6  8]

V[3]:=17:

m1[1,2]; # V and m1 do not share data

                                      3

m2[1,2]; # the difference, V and m2 share data

                                     17

whattype(m1), whattype(m2);

                               Matrix, Matrix

acer

Using the GlobalOptimization package (under 64bit Maple 14.01 on Windows 7),

restart:
with(GlobalOptimization):

GlobalSolve(
   (3*a/sqrt(1-u^2)+b/sqrt(1-t^2))/c,
   {a^2+2*b*c*u >= b^2+c^2, a*t+b*u <= c,
   b^2*(t^2-u^2)/(t^2-1)+c^2 <= 2*b*c*u},
   t=0..1, u=0..1, a=0..100, b=0..100, c=0..100);

   [4.00000000000000711, [a = 35.3564579217052, b = 35.3563778458042, 

     c = 50.0015242889233, t = 0.707107562133025, u = 0.707105960663738]]

That was for the default 'multistart' method. Unfortunately it's not clear how it chooses initial points. As the following userinfo messages indicates it appears to have already obtained a good result prior to the userinfo printouts. These messages seem to indicate that it has found the objective value very close to 4.0000 during some initial local search phase, before any global search phase has begun. Notice that the default feasibility tolerance used by GlobalSolve is 1.0e-4 which is really quite coarse a possible violation. But I was also able to get the same result from `GlobalSolve` using feasibility=1e-11, but no finer. The regular `Optimization:-Minimize` command can also find the same objective value if supplied with that same coarse feasibility tolerance (but doesn't get nearly as close to 4.00000000000000 for feasibilitytolerance=1e-5 unless the iterationlimit is supplied as high as 400).

restart:
with(GlobalOptimization):
infolevel[GlobalOptimization]:=100:
GlobalSolve(
   (3*a/sqrt(1-u^2)+b/sqrt(1-t^2))/c,
   {a^2+2*b*c*u >= b^2+c^2, a*t+b*u <= c,
   b^2*(t^2-u^2)/(t^2-1)+c^2 <= 2*b*c*u},
   t=0..1, u=0..1, a=0..100, b=0..100, c=0..100,
   method=singlestart);
GlobalSolve: calling NLP solver
SolveGeneral: calling global optimization solver
SolveGeneral: number of problem variables 5
SolveGeneral: number of nonlinear inequality constraints 3
SolveGeneral: number of nonlinear equality constraints 0
SolveGeneral: method singlestart
SolveGeneral: merit function evaluation limit 8000
SolveGeneral: non-improving merit function evaluation limit 1600
SolveGeneral: constraint penalty multiplier 100.0
SolveGeneral: target merit function value -0.10e11
SolveGeneral: local search target objective function value -0.10e11
SolveGeneral: local search feasibility tolerance 0.10e-5
SolveGeneral: local search optimality tolerance 0.10e-5
SolveGeneral: time limit in seconds 800
SolveGeneral: trying evalhf mode
Extcall: Lower bounds, nominal values, and upper bounds of decision variables:
Extcall: 1   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 2   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 3   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 4   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: 5   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: --- Initial function evaluation at nominal solution ---
Extcall: Objective function value 4.61880215351700585
Extcall: 1   5.000000e+001
Extcall: 2   5.000000e+001
Extcall: 3   5.000000e+001
Extcall: 4   5.000000e-001
Extcall: 5   5.000000e-001
Extcall: Constraint function values, 1, 0.
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, 0.
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Number of function evaluations: 259
Extcall:  Current optimum estimate: 4.00000000000000533
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: 1   -9.094947e-013
Extcall: 2   0.000000e+000
Extcall: 3   4.547474e-013
Extcall: --- Automatic Global Optimization Procedure ---
Extcall: --- Global Adaptive Random Search ---
SolveGeneral: trying evalf mode
Extcall: Lower bounds, nominal values, and upper bounds of decision variables:
Extcall: 1   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 2   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 3   0.000000e+000  5.000000e+001  1.000000e+002
Extcall: 4   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: 5   0.000000e+000  5.000000e-001  1.000000e+000
Extcall: --- Initial function evaluation at nominal solution ---
Extcall: Objective function value 4.61880215351700585
Extcall: 1   5.000000e+001
Extcall: 2   5.000000e+001
Extcall: 3   5.000000e+001
Extcall: 4   5.000000e-001
Extcall: 5   5.000000e-001
Extcall: Constraint function values, 1, 0.
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, 0.
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Number of function evaluations: 259
Extcall:  Current optimum estimate: 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: 1   -1.364242e-012
Extcall: 2   0.000000e+000
Extcall: 3   1.818989e-012
Extcall: --- Automatic Global Optimization Procedure ---
Extcall: --- Global Adaptive Random Search ---
Extcall: The relative objective function improvement isless than .100000000000000002e-7 in 3 subsequent major iterations.
Extcall: The local search phase is terminated
Extcall: Best solution found in GARS + LS search mode
Extcall: Number of function evaluations: 8398
Extcall:  Current optimum estimate: 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: Constraint function values, 1, -.136424205265939236e-11
Extcall: Constraint function values, 2, 0.
Extcall: Constraint function values, 3, .181898940354585648e-11
Extcall: --- Summary of Optimization Results ---
Extcall: Total number of function evaluations: 8399
Extcall: Estimated optimum (including penalties) 4.00000000000000711
Extcall: Estimated optimal solution vector (components):
Extcall: 1   3.535646e+001
Extcall: 2   3.535638e+001
Extcall: 3   5.000152e+001
Extcall: 4   7.071076e-001
Extcall: 5   7.071060e-001
Extcall: All stated constraints are met within tolerance .99999999999999995e-6
Extcall: Sum of constraint violations: .181898940354585648e-11
Extcall: Solution time (seconds), including user and system I/O operations, 0.
SolveGeneral: total number of function evaluations 8399
SolveGeneral: runtime in external solver 0.
SolveGeneral: maximum constraint infeasibility 0.181898940354585648e-11
SolveGeneral: cycling or stall detected in solver

       [4.00000000000000711, [a = 35.3564579217052, b = 35.3563778458042, 
                              c = 50.0015242889233, t = 0.707107562133025,
                              u = 0.707105960663738]]

acer

Does it help at all, to set u=t=sqrt(2)/2 ?  Just curious.

What about if you search on a=0.0..0.0001, b=0.0..0.0001, c=0.0..0.001 ?

acer

@Subjectivity Here's an attempt.

Download noiseweight.mw

What hapens if you "declare" the parameters of procedure `iang1`, using double-colon?

I mean, something like,

iang1 := proc (m1::float[8], vF1::float[8], omega::float[8], q::float[8]) ....

This would help the Compiler know that these were not going to be integer, with whatever other effects would ensure for that in the code (casts, etc).

acer

You want S.E.T to be approximately diagonal, right? (...by defn of Smith Form)

If so, then you'll need to increase the working precision (ie. `Digits`). But as the working precision gets higher, the rank of floating-point Matrix E may change. (The small entries of E can become significant, as the working precision increases.)

restart:
interface(rtablesize=100):

E := ImportMatrix(cat(kernelopts(homedir),
                      "/Downloads/model15_E.txt"),
                  delimiter=" "):

with(LinearAlgebra):

for i from 10 to 17 by 1 do
   Digits:=i;
   S,SNF,T := evalf(SmithForm(E,output=['U','S','V']));
   printf("   %ld      %3.3g     %.3e   %ld   %ld\n",
          Digits, Determinant(E), Norm(S.E.T - SNF),
          Rank(SNF), Rank(E));
end do:

for i from 10 to 17 do
   Digits:=i;
   S,SNF,T := evalf(SmithForm(E,output=['U','S','V'])):
   Z:=S.E.T:
   printf("Digits=%ld\n",i);
   printf("%.3f\n\n\n",map(fnormal,Z));
end do:

acer

You've defined `r` as an expression, but you use it as if it were a parameter name in your call to the plot3d command.

After assigning to your constants, the indeterminate names left in `Top` and `Bottom` are `phi` and `v`.

Did you instead mean something like,

plot3d([Top, Bottom], phi = s .. 20, v = 0 .. 2*Pi);

What coordinate system are you trying to plot in?

What are these two lines supposed to accomplish? (nb. `u` doesn't even appear anywhere else.)

u := u1+(u2-u1)*JacobiSN((1/2)*phi*sqrt(s*(u3-u1)), k)^2; 

r := 1/u;

acer

Try it with `algsubs` instead of `subs`.

algsubs(y2*y3 = 0, u(y2, y3));

                             2     2
                           y2  + y3 

acer