Carl Love

Carl Love

26104 Reputation

25 Badges

11 years, 57 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

This is a good question. I'm suprised noone has answered it yet.

There is no need to take the derivative inside the procedure. Symbolically, the steps for doing so are the same for any j. By putting it in a procedure, you are forcing it to go through the steps of taking essentially the same derivative every time the procedure is called.

G := (1/2)*(YY[j]-YY[j-1])^2/(XX[j]-XX[j-1])+(1/2)*(YY[j+1]-YY[j])^2/(XX[j+1]-XX[j]);

Gp:= diff(G, XX[j]);

(1/2)*(YY[j]-YY[j-1])^2/(XX[j]-XX[j-1])+(1/2)*(YY[j+1]-YY[j])^2/(XX[j+1]-XX[j])

-(1/2)*(YY[j]-YY[j-1])^2/(XX[j]-XX[j-1])^2+(1/2)*(YY[j+1]-YY[j])^2/(XX[j+1]-XX[j])^2

Now here I evaluate this for specific lists XX and YY, and specific j. I used lists of symbols so that you can see what's going on. But you can use lists (or vectors) of numbers.

X:= [a,b,c,d]:  Y:= [w,x,y,z]:

eval(Gp, [XX= X, YY= Y, j= 2]);

-(1/2)*(x-w)^2/(b-a)^2+(1/2)*(y-x)^2/(c-b)^2

If you want to make this into a procedure, now is the time, and do it with unapply.

GX:= unapply(Gp, [XX,YY,j]);

proc (XX, YY, j) options operator, arrow; -(1/2)*(YY[j]-YY[j-1])^2/(XX[j]-XX[j-1])^2+(1/2)*(YY[j+1]-YY[j])^2/(XX[j+1]-XX[j])^2 end proc

GX(X,Y,2);

-(1/2)*(x-w)^2/(b-a)^2+(1/2)*(y-x)^2/(c-b)^2

 

 

Download FinDiff.mw

By transforming the problem to spherical coordinates, we can reduce the number of independent variables from three to two; this also makes plotting the milieu much easier. A well-constructed plot will lead us directly to the absolute minimum without any ambiguity about local minima.



 

restart;

Digits:= 15:

The objective function:

A:= x^2*y^2/2 + y^2*z^2 + z^2*x^2 + 96/(x+y+z+1);

(1/2)*x^2*y^2+y^2*z^2+z^2*x^2+96/(x+y+z+1)

Note that this is symmetric in x and y.

 

The constraints (or domain) is the sphere of radius sqrt(5) centered at the origin, in the first octant. Knowing that, the constraint equation is not used for the ret of this problem.

 

Set up spherical coordinate transformation.

r:= sqrt(5):

Trans:= [x,y,z] =~ r*~[(cos(theta),sin(theta))*~sin(phi), cos(phi)];

[x = 5^(1/2)*cos(theta)*sin(phi), y = 5^(1/2)*sin(theta)*sin(phi), z = 5^(1/2)*cos(phi)]

Convert the constrained objective to spherical coordinates.

AS:= eval(A, Trans);

(25/2)*cos(theta)^2*sin(phi)^4*sin(theta)^2+25*sin(theta)^2*sin(phi)^2*cos(phi)^2+25*cos(phi)^2*cos(theta)^2*sin(phi)^2+96/(5^(1/2)*cos(theta)*sin(phi)+5^(1/2)*sin(theta)*sin(phi)+5^(1/2)*cos(phi)+1)

Now we have only two independent variables to deal with, phi and theta, rather than three: x, y, z.

 

Plot the domain using the objective function to color the surface.

plot3d(
   [r,theta,phi], theta= 0..Pi/2, phi= 0..Pi/2,
   coords= spherical,
   color= AS,
   orientation= [45,45,0], grid= [25,25], thickness= 0,
   axes= boxed
);

(This plot is much clearer in actual than it appears on MaplePrimes.) The values of the objective function, low to high, correspond to the color spectrum, red to violet. Knowing that, it is clear that there are three local minima---the two red spots and the one yellow spot---and that one of the red spots is the absolute minumum (red being less than yellow). Because of the x-y symmetry, the objective value at each red spot is the same. The situation for maxima is not so clear: They all lie on the border. But we aren't interested in maxima. By counting the grid lines, it is clear that one of the red spots occurs close to theta = 1/4*Pi/2, phi = 22/25*Pi/2.

 

Minimize the objective by solving for where its gradient is 0.

Sol:= fsolve({diff(AS,phi),diff(AS,theta)}, {theta=1/4*Pi/2, phi= 22/25*Pi/2});

{phi = 1.36079017176763, theta = .414107150255695}

And the minimal value is

evalf(eval(AS, Sol));

24.6692515212380

At (x, y, z) =

evalf(eval(Trans, Sol));   

[x = 2.00209163810771, y = .879965268827956, z = .466143967327395]

 

 



Download Min_on_sphere.mw

 

A shortcut for the transpose:

B^%T;

This works even if LinearAlgebra has not been loaded.

(1) The csgn is just a function that is 1 or -1 depending on the sign of the (real part of) its argument. (It's slightly more complicated than that. See ?csgn for full details.) You can eliminate its appearance in your answer by making some assumptions. I made the assumptions R>0, H>z, z>0, and I got a much simplified answer for the integral; but I don't know if these assumptions are appropriate for your problem.

(int(int(theta*R*(H-p)/(H*(L*L)*sqrt((R*(H-p)/H)^2+(z-p)^2)), p = 0 .. z), o = 0 .. 2*Pi))/(4*Pi*epsilon) assuming R>0, H>z, z>0;

(2) It's better to use unapply to make V into a procedure than to put the integral on the right side of -> (the operator-defining arrow):

J:= int(theta*R*(H-p)/(H*(L*L)*sqrt((R*(H-p)/H)^2+(z-p)^2)), p = 0 .. z);

V:= unapply(1/2/epsilon*J, z);

If you keep it the way you had it, then the integral gets redone every time that fsolve evaluates V. By using unapply, the integral is only done once.

(3) You must have supplied numeric values for theta, H,and L that you don't show,because fsolve would have returned an error otherwise. Anyway, the appearance of csgn is nothing to worry about. It's only a 1 or -1.

 

 

All levels of algebraic grouping in Maple are done with ordinary parentheses (), not square brackets [] or curly braces {}.

Your file has an extra blank line at the end. Use an ordinary text editor, such as Windows Notepad, to delete that. After doing that, Maple's Import Data Assistant works perfectly. Enter

theData:= ImportData();

In the ensuing dialog, select your file, and for the file type, select "Delimited". You can accept the defaults for the rest of the dialog, which guesses correctly that the delimiter is ;, that the first 6 lines should be ignored, etc.

It's surprising that the Assistant is so good at trimming extra junk from the beginning of the file, but is thrown off by a simple extra blank line at the end.

If you actually want to read in the string column headers, please respond, and we can do that separately. Surely you want the numerical data separated as a Matrix.

I don't think that any of the respondents here has yet pointed out explicitly that the original input exponent, 29.403243784, has 11 significant digits. Maple's default working precision is 10 digits, so the first internal step in the computation is to round this to 29.40324378. Using this as the exponent, Maple's result is accurate to the full 10 digits, well within the required 0.6 ulps. This is why the major improvement in accuracy occurs when going from 10 to 11 Digits, as acer showed.

So, to answer explicitly your question: To avoid this problem, set the global variable Digits to a value at least as high as the number of significant digits in the most precisely specified operand. Note that significant digits is not the same thing as digits after the decimal point.

Digits:= 11;

You did not specify an initial value for n. I will assume it's 0. Your function u is called Heaviside in Maple (and elsewhere).

y:= t-> Sum(sin(t-n*Pi)*Heaviside(t-n*Pi), n= 0..infinity);

plot(y, 0..30);

Here is another solution. It is similar to Kitonum's in that it iterates over all possible a, b, c; but it stops iterating the recursive sequence as soon as a noninteger is found. This is done by checking the remainder of the division rather than with a type check. The solutions are stored in a table rather than by the asymptotically quadratic list-appending technique (hardly makes any difference here since there are only two solutions, but it is good to practice the avoidance of list appending). All in all, there's about a factor-of-5 performance improvement, mostly from the early terminations of the sequence iteration.

restart;
t:= Array(-1..20, [0,1]):  # [0,1] are initial values
Soln:= table():
abc:= combinat:-cartprod([[$1..20] $ 3]):  #Iterator
while not abc[finished] do
    assign(('a','b','c')= abc[nextvalue]()[]);
    for n from 0 to 19 do
        t[n+1]:= iquo((a*n*(n+1)+b)*t[n] + c*n^2*t[n-1], (n+1)^2, 'r');
        if r <> 0 then  break  end if
    end do;    
    if r = 0 then  Soln[a,b,c]:= convert(t,list)  end if
end do:
eval(Soln);

You have four specific x values in you IBCs, but Maple can handle at most two. You have -5000, -1000, 1000, 5000.

I don't know why the surd doesn't work, but limit(a^(1/k), k= infinity) does work.

There seems to be a bug when there's only one element in the list being mapped over. The bug is not about a difference between `^` and op; but rather about a difference between [1,2] and [[a,b]]: The former has two elements while the latter has only one. The bug is also manifested by Map(`^`, [2], 2).

Here's a very basic (although structured, functional, and extensible) implementation. I didn't include any error checking of user input. So the user needs to make sure that there's the right number (25 by default) of characters in the key, split doubletons in the plain text, make sure only alphabet characters are in the plaintext or ciphertext, and make sure that texts have an even number of characters. However, I did generalize it so that you are not restricted to (the traditional) 5x5 alphabet: You can use any rectangle (just change R and C in the module locals).

(*
A basic implementation of the Playfair cipher (see "Playfair cipher" in Wikipedia).
There is no error checking of user input, and the user has to split doubletons on
their own.
                                                                                   *)
Playfair:= module()
option `Written 19-Apr-2013 by Carl Love.`;
local
    R:= 5, C:= 5,  #For a 5x5=25 alphabet
    
    #Convert number in 1..25 into row-column coordinates
    to_row_col:= proc(n)  local r,c;  r:= iquo(n-1, C, 'c');  (r+1, c+1)  end proc,

    #Convert row-column coordinates to number in 1..25
    from_row_col:= (r,c)-> C*(r-1)+c,

    rotf:= k-> k mod C + 1,  #mod C, shifted by 1 (rotate forward)
    rotd:= k-> k mod R + 1,  #mod R, shifted by 1 (rotate down)
    
    #For a pair in 1..R*C, return the corresponding Playfair pair. Note that if
    #we just consider the indices in the array, then the mapping from index-pair
    #to index-pair is a constant, independent of the cipher key.
    Playfair_pair:= proc(a,b)
    option remember;
    local r1, r2, c1, c2;
        if a=b then  error "expected unequal arguments"  end if;
        (r1,c1):= to_row_col(a);
        (r2,c2):= to_row_col(b);
        if r1=r2 then  (from_row_col(r1,rotf(c1)), from_row_col(r2,rotf(c2)))
        elif c1=c2 then  (from_row_col(rotd(r1),c1), from_row_col(rotd(r2),c2))
        else (from_row_col(r1,c2), from_row_col(r2,c1))
        end if
    end proc,

    EnDig:= table(),  #Encrypt: Digraph to Digraph
    DeDig:= table()   #Decrypt: Digraph to Digraph
;
export
    SetKey:= proc(Key::string)
    local i, j, i2, j2, DigPlain, DigCrypt;
         #Encrypt and decrypt every possible digraph.
        for i to R*C do  for j to R*C do
            if i=j then  next  end if;
            (i2,j2):= Playfair_pair(i,j);
            DigPlain:= cat(Key[i],Key[j]);
            DigCrypt:= cat(Key[i2],Key[j2]);
            EnDig[DigPlain]:= DigCrypt;
            DeDig[DigCrypt]:= DigPlain
        end do  end do;
        [][]
    end proc,

    Encrypt:= proc(plaintext::string)
    local k;
        cat(seq(EnDig[plaintext[2*k-1..2*k]], k= 1..iquo(length(plaintext),2)))
    end proc,
    
    Decrypt:= proc(ciphertext::string)
    local k;
       cat(seq(DeDig[ciphertext[2*k-1..2*k]], k= 1..iquo(length(ciphertext),2)))
    end proc
;
end module;

 

restart;

 

(*

module () local R, C, to_row_col, from_row_col, rotf, rotd, Playfair_pair, EnDig, DeDig; export SetKey, Encrypt, Decrypt; end module

Here is the example from the Wikipedia page "Playfair cipher"

Playfair:-SetKey("PLAYFIREXMBCDGHKNOQSTUVWZ"); #No J

Note that X is inserted in "tree" to split the double e.

Playfair:-Encrypt("HIDETHEGOLDINTHETREXESTUMP");

Playfair:-Decrypt(%);

 

 

Download Playfair.mw

The package linalg, the command evalm, and the operators &. and &* (when used in the linalg/matrix context) should not be used... ever. These commands are ancient and deprecated. Even experts are forgetting the nuances of their usage. If you have a professor giving you code to work with, please tell that person that I said to not use these commands... ever.

Maple's solve does not return anything when there are no solutions. I agree that that is confusing. In this case, some trivial algebra "by hand" shows that there is indeed no solution. The -2s and rs cancel and you're left with

abs(1+sqrt(1+4*r)) < 1

On the left, you have 1 plus something nonnegative, so the inequality is never true.

The assuming has nothing to do with this: There are no solutions no matter what you assume.

First 358 359 360 361 362 363 364 Last Page 360 of 378