Carl Love

Carl Love

28020 Reputation

25 Badges

12 years, 302 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Here's how to do it with the order. That is, how to run Maximize and save each call to the objective function and the order that the calls were made.

 

(**)

restart:

Here's an example maximization problem.

(**)

Optimization:-Maximize(sin(x), {x >= 0, x <= 3*Pi/2});

[1., [x = HFloat(1.5707963267950427)]]

Now I will convert this problem into a form where we save each call to the objective function for later use.

(**)

Problem:= module()
local
     ct,
     Calls
;
export
     Obj:= proc(x)
     local r:= sin(x);  
          ct:= ct+1;  
          Calls[ct]:= (x,r);  
          r  
     end proc,

     Constr1:= x-> 0 - x,  # x >= 0
     Constr2:= x-> x - 3*Pi/2,  # x <= 3*Pi/2
     Constrs:= {Constr1, Constr2},

     Init:= proc()  
          ct:= 0;  
          Calls:= table()
     end proc,

     AllObjEvals:= proc()  
     local k;  
          [seq([Calls[k]], k= 1..ct)]
     end proc
;
     Init()
end module:

 

(**)

Optimization:-Maximize(Problem:-Obj, Problem:-Constrs);

[1., Vector(1, {(1) = 1.57079632679490})]

(**)

Data:= Problem:-AllObjEvals();

[[HFloat(0.0), HFloat(0.0)], [0., 0.], [0.50000000000000000000000000e-8, 0.49999999999999999791666667e-8], [-0.50000000000000000000000000e-8, -0.49999999999999999791666667e-8], [HFloat(1.0), HFloat(0.8414709848078965)], [1., .84147098480789650665250232], [1.0000000050000000000000000, .84147098750940802547481359], [.99999999500000000000000000, .84147098210638496679341644], [HFloat(2.1753426496700228), HFloat(0.8227600473221269)], [2.175342649670023, .82276004732212685306409917], [2.1753426546700230000000000, .82276004448018250601607050], [2.1753426446700230000000000, .82276005016407117954312665], [HFloat(1.569745684590937), HFloat(0.9999994480755304)], [1.569745684590937, .99999944807553039954776143], [1.5697456895909370000000000, .99999944808078359710110782], [1.5697456795909370000000000, .99999944807027717699442884], [HFloat(1.570855739119737), HFloat(0.9999999982350878)], [1.570855739119737, .99999999823508782904968963], [1.5708557441197370000000000, .99999999823479075492566251], [1.5708557341197370000000000, .99999999823538487817371679], [HFloat(1.5707963267845844), HFloat(1.0)], [1.570796326784584, .99999999999999999999994682], [1.5707963317845840000000000, .99999999999999998755150992], [1.5707963217845840000000000, .99999999999999998744838373], [HFloat(1.570796326794897), HFloat(1.0)], [1.570796326794897, 1.0000000000000000000000000], [1.5707963317948970000000000, .99999999999999998749999810], [1.5707963217948970000000000, .99999999999999998750000190]]

(**)

 

Please let me know if you are able to apply this to your problem.

Download Max_with_saved_dat.mw

conjugate(1/sqrt(I*t*lambda+m))*sqrt(lambda^2*(t^2)+m^2)*(1/sqrt(I*t*lambda+m)):
simplify(evalc(%));
                               1

It would not be possible for Maple in general to implement what you want. Maple deals with polynomial arithmetic over finite fields. x^2+x+1 is considered an irreducible polynomial over Z[2]; it is not equivalent to 1. So Z[2]/(x^2+x+1) is the finite field with 4 elements.

But, you can easily redefine mod to suit your purposes:

restart:
Modp:= eval(modp):
unprotect(modp):
modp:= (e,m)-> Modp(`if`(m=2, evalindets(e, `^`, x-> op(1,x)), e), m):
protect(modp):

x^2 mod 2;
                                            x

a*b*c+a^2*c*d+d+c^2*d^2*a^2*b^2+d^2*a*b^2*c+d^2*a*b^2+a*c+c*d+d*b+d*a^2*b*c+c+c^2+d^2*a^2*c+d^2*a^2*b^2*c+d^2*b^2+d^2*a+
d^2+b*c^2+c^2*d^2+b*c^2*d+a^2*c^2*b*d+c^2*d^2*a^2+c*d*a*b+a*b^2*c*d+b^2*c*d+b^2*c^2*d+d^2*a*c+a^2*b^2*c^2*d+b*c+d*a*b+d*a+
a^2*c^2+a^2*c^2*b+c^2*d^2*b^2 mod 2;
                                             0

Yes, function application distributes over most other operators; and constants, when applied as functions, become constant functions. Example:

((D-2)^3)(sin)(x);

There is a missing character between your k and 1. I will assume it is a plus sign. You can enter the expression pretty much as you currently have it into Maple:

Limit((Sum((n-k)/(k+1), k= 0..n-1) -ln(n!))/n, n= infinity);

value(%);

                                      γ

The answer is hard to read when I cut-and-paste it above, but it is the Euler-Mascheroni constant, gamma.

 

No, it can't be done.

Please post the example. It might help if you give fsolve a hint by making an initial guess based on an implicitplot.

See the section "Fundamental Recurrence Relation" of the Wikipedia article "Generalized Continued Fraction".


(**)

restart:

(**)

recA:= A(n+1)=b(n+1)*A(n) + a(n+1)*A(n-1), A(-1)=1, A(0)=b(0):

(**)

recB:= B(n+1)= b(n+1)*B(n) + a(n+1)*B(n-1), B(-1)=0, B(0)=1:

(**)

b:= n-> 2:

(**)

a:= n-> (2*n-1)^2:

(**)

An:= rsolve({recA}, A(n));

2^(n+1)*GAMMA(n+3/2)*(1/Pi^(1/2)-(1/4)*((-1)^(n+1)*Psi(5/4+(1/2)*n)-(-1)^(n+1)*Psi(3/4+(1/2)*n)-Pi)/Pi^(1/2))

(**)

Bn:= rsolve({recB}, B(n));

-(1/4)*2^(n+1)*GAMMA(n+3/2)*((-1)^(n+1)*Psi(5/4+(1/2)*n)-(-1)^(n+1)*Psi(3/4+(1/2)*n)-Pi)/Pi^(1/2)

(**)

simplify(An/Bn) assuming n::posint;

(2*(-1)^n*Psi((1/2)*n+1/4)*n+2*(-1)^(n+1)*Psi(3/4+(1/2)*n)*n+2*Pi*n+(-1)^n*Psi((1/2)*n+1/4)-(-1)^n*Psi(3/4+(1/2)*n)+Pi+4*(-1)^n+8*n+4)/(2*(-1)^n*Psi((1/2)*n+1/4)*n+2*(-1)^(n+1)*Psi(3/4+(1/2)*n)*n+2*Pi*n+(-1)^n*Psi((1/2)*n+1/4)-(-1)^n*Psi(3/4+(1/2)*n)+Pi+4*(-1)^n)

(**)

limit(%, n= infinity);

(4+Pi)/Pi

(**)

simplify(4/(1+%));

2*Pi/(Pi+2)

(**)

 


Download gen_cont_frac.mw

It's as simple as

(A_old,A_new):= (A_new,A_old):

Ordinary assignment (via :=) only copies the addresses of tables and arrays, which is what you want. Multiple ordinary assignment is done in such a way that an intermediate variable is not necessary.

1. You probably used an e somewhere earlier in your computation where you had meant to use exp. If you print the expression with lprint, it will show you if the e is actually exp(1).

2. Yes, the argument of a positive real number is 0.

It is surprising that there is no single command to do this. So I wrote a procedure to do it.

To use this, you'll need to download a combinatorics package called Iterator from the Maple Applications Center. The package has a command SetPartitions which lets you control the size of the blocks in the partition, but you can't directly control the number of blocks. So I iterate over the integer partitions of n to get all possible combinations of block sizes with k blocks. Then these need some modification to be put into the form expected by SetPartitions.

PartitionIntoKBlocks:= proc(n::{set,nonnegint}, k::nonnegint)
option `Written by Carl Love 2013-Sep-27.`;
uses J=Iterator, C= combinat, S= Statistics;
local s, P, N:= `if`(n::set, n, {$1..n});
     {seq(
          seq({s[]}, s= J:-SetPartitions(N, map(`[]`@op, S:-Tally(P)), compile= false)),
          P= select(x-> nops(x)=k, C:-partition(nops(N)))
     )}
end proc:

PartitionIntoKBlocks(4,2);
    {{{1}, {2, 3, 4}}, {{2}, {1, 3, 4}}, {{3}, {1, 2, 4}},
      {{4}, {1, 2, 3}}, {{1, 2}, {3, 4}}, {{1, 3}, {2, 4}},
      {{1, 4}, {2, 3}}}

Joe Riel, if you're reading: I think Iterator:-SetPartitions should have an option for the number of blocks in the partition.

This was trickier than I expected, probably because of a bug in Maple.

a,b,c:= 'RootOf(x^3-x^2-x-1, x, index= k)' $ k= 1..3:
evala(
     evala((a^2014-b^2014)/(a-b)) +
     evala((b^2014-c^2014)/(b-c)) +
     evala((c^2014-a^2014)/(c-a))
);

8334826359539004144901989461804527198477505526542440583547585702\
  52333926756888184793278322997734194204312479082099178141286346\
  19093619632692976466847552803713809622682218736343227666160495\
  08221814701281712551536720818416541447585602970496630535429027\
  14035268777406803186312918798664076365740099763393011180604524\
  46886004153088055743393327297340352231807305287646640702341995\
  05000473477741234080182885893584028536441114628416893583917671\
  56862792958972093429609456884388481631837778020024880944176069\
  54812869711446457823816195550797524

The answer is just a large integer. It is a[2014] for the Fibonacci-like recurrence a[n] = a[n-1]+a[n-2]+a[n-3], a[0]=0, a[1]=3, a[2]=2.

The bug ensues if you try to do it with a single evala command. It seems to get stuck in an infinite loop.


(**)

restart:

Note that Maple's series command will return the linear and quadratic approximations. (The taylor command does the same thing.)

Linear is considered order 2:

(**)

convert(series(f(x), x=a, 2), polynom);

f(a)+(D(f))(a)*(x-a)

Quadratic is order 3:

(**)

convert(series(f(x), x=a, 3), polynom);

f(a)+(D(f))(a)*(x-a)+(1/2)*((D@@2)(f))(a)*(x-a)^2

Applying this to your function, we get:

(**)

P1:= convert(series(arccos(x), x=0, 2), polynom);

(1/2)*Pi-x

(**)

P2:= convert(series(arccos(x), x=0, 3), polynom);

(1/2)*Pi-x

They are the same because x=0 is an inflection point, which will be clear from the plot.

(**)

plot([arccos(x), P1, P2], x= -1..1);

(**)

 


Download linear_approx.mw


(**)

restart:

(**)

Digits:= 2^6:

(**)

eq:= 6.750969101*10^26*_Z*1000^(9/59)*800^(41/59)-1.015701668*10^27*1000^(50/59)-1.189275761*10^27*800^(50/59)+1.103397224*10^26*_Z*1000^(32/59)*800^(18/59)-7.457053114*10^23*3^(1/5)*1000^(41/59)*800^(50/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-7.457053114*10^23*3^(1/5)*1000^(50/59)*800^(41/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-1.988547497*10^26*3^(1/5)*1000^(32/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)-2.485684371*10^26*3^(1/5)*800^(32/59)*_Z^(4/5)*sqrt(12500)*4^(19/20)*(800^(9/59)+1000^(9/59))^(1/50)*(702.222720*1000^(9/59)+976.508160*800^(9/59))^(4/5)+5.575397080*10^26*_Z*1000^(41/59)*800^(9/59)+1.308382216*10^26*_Z*1000^(18/59)*800^(32/59)+1.378754484*10^27*_Z*800^(50/59)+1.114796623*10^27*_Z*1000^(50/59)-1.144263788*10^26*1000^(32/59)*800^(18/59)-1.028573111*10^26*1000^(18/59)*800^(32/59)-5.399913635*10^26*1000^(41/59)*800^(9/59)-5.544622178*10^26*1000^(9/59)*800^(41/59):

(**)

indets(eq, {name, name^numeric});

{_Z, _Z^(4/5)}

Use a change of variables to convert to a polynomial:

(**)

eq_simp:= simplify(eval(eq, _Z= z^5), symbolic);

1252302292256880108322573586627.994134196842318012573369627002321*z^5-1113835100487674276481580185899.975871697792466972768109672627414-33206779716982073462031352179419335.13149299607560858653432559958*z^4

(**)

r:= fsolve(eq_simp);

26516.58463160465918042279703097330382230532550324190250597416872

Residual in the simplified equation:

(**)

eval(eq_simp, z= r);

0.

Resisual in the original equation:

(**)

evalf(eval(eq, _Z= r^5));

-0.324e-10

(**)

sol1:= r^5;

13109554349246137246800.44849539580049816642933649202738691745839

For comparison, apply RootFinding:-Analytic to the original equation, using the same bounding box as did Markiyan. (This takes much longer than fsolve.)

(**)

sol2:= RootFinding:-Analytic(eq, _Z= 0-10*I..10^30+100*I);

13109554349246137246800.44849539580049816642933649202738691745852

(**)

sol1-sol2;

-0.13e-39

A plot confirms this root.

(**)

plot(eq, _Z= 0..2*10^22);

(**)

 


Download sensitive_equatio.mw

First, it helps tremendously if you indent your code, like this:

restart:
P := array([[1, 4], [2, 3], [3, 2], [4, 1], [6, 5], [6, 1], [6, 2], [5, 3]]);

j := 1; k := j+2; Fold := 0; Fnew := 1; counter := 0;

for i to 6 do
     j := i; k := j+2;
     while Fnew > Fold do
          Fold := Fnew;
          L[j, k] := evalf((P[k, 1]-P[j, 1])^2+(P[k, 2]-P[j, 2])^2)^(1/2);
          m[j, k] := (P[k, 1]-P[j, 1])/(P[k, 2]-P[j, 2]);
          a := m[j, k];
          b := -1;
          c := P[j, 2]-m[j, k]*P[j, 1];
          for i from j+1 to k-1 do
               e[i] := evalf(abs(a*P[i, 1]+b*P[i, 2]+c)/sqrt(a^2+b^2));
               E[j, k] := sum(e[m], m = j+1 .. k-1)
          end do;
          Fnew := L[j, k]-E[j, k];
          counter := counter +1;
          k := k+1
     end do;
     right[j] := k-2
 end do;

And when I look at it like that, I immediately see the problem: You used the same index variable, i, for both the inner and outer for loops. So I simply changed the loop index to ii for the inner loop, and I changed its three invocations within the loop to ii also:

          for ii from j+1 to k-1 do 
               e[ii] := evalf(abs(a*P[ii, 1]+b*P[ii, 2]+c)/sqrt(a^2+b^2));
               E[j, k] := sum(e[m], m = j+1 .. k-1)
          end do;

And then I got 6 values in right.

First 345 346 347 348 349 350 351 Last Page 347 of 395