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

When you write of testing inequality, do you mean that you want it to return 'false' if at least one pair of entries violates the condition? (That's as opposed to all pairs having to violate it.)

If so then it can be faster to have a process which quits upon the first elementwise violation and which does not necessarily have to walk all entries. (I'm sure Joe Riel can beat whatever I come up with, soundly, for speed.)

I don't think that there are commands `orzip` and `andzip`.

Suppose that we wish to test whether A ≥ B, and by which we have taken the meaning to be that this is true only if A[i,j]≥B[i,j] for all i,j.

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   if ormap(`<`,A-B,0) then
      # restriction was false, for at least one pair of entries
      count:=count+1;
  end if;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 8.767 seconds

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   try
      zip(proc(a,b) if a < b then error "found"; end if; end proc, A, B);
   catch "found":
      # restriction was false, for at least one pair of entries
      count:=count+1;
   end try;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 0.671 seconds

restart:
randomize(4567890):
with(LinearAlgebra):
A:=RandomMatrix(4,4,generator=10..100):
st:=time():
count:=0:
for i from 1 to 10000 do
   B:=RandomMatrix(4,4,generator=-80..20);
   if ormap(evalb,`<`~(A,B)) then
      # restriction was false, for at least one pair of entrie
      count:=count+1;
  end if;
end do:
printf("\n     %a iterations produced false in %.3f seconds\n\n",count,time()-st);

     887 iterations produced false in 0.405 seconds

acer