Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Simply,

Gcd(f,g) mod 3;

Note that the capital G distiguishes this from the "regular" command gcd. See the help at ?mod. This will have a list of all the commands designed to work with mod over finite rings.

What's immediately obvious to me is that you're not updating X and Y. I think that you think that your updates of R implicitly also update X and Y, but that doesn't happen. They can only be updated when they appear on the left side of an assignment statement.

My advice is to not use R at all; just maintain X and Y. That is, for each of the three values of r there should be two assignment statements: one for X and one for Y.

This is one of the most common Maple errors, called premature evaluation. The functions I1 and I2 obviously can't be evaluated numerically until t has a numeric value. In the plot command you invoke I1(t) and I2(t). These invocations are evaluated before t is given any numeric value. One solution is this:

plot([I1,I2], 0..200, color= [red, blue]);

P.S.: The above is a correction to the premature evaluation problem only. I haven't tested whether this particular plot actually works. Markiyan got a nonreal value. If you get too many nonreal values, the plot command will give you an empty plot.

Here is an example random problem with 15 variables, nine general constraints, and all variables restricted to the unit interval:

V:= [x||(1..15)]:
Objective:= randpoly(V, degree= 1, dense):
Constraints:= ['randpoly(V, degree= 1) >= 0' $ 9]:
Optimization:-LPSolve(Objective, Constraints, (V =~ 0..1)[]);

This code executes instantaneously. However, you specified the open interval 0 < x < 1. Is that what you really want, as opposed to the closed interval 0 <= x <= 1? Open-interval constraints can't be done by the program, and it isn't even clear whether they can be solved at all. Continuous functions are guaranteed to have extrema on closed intervals, not on open intervals. You can get an approximate solution by choosing a small positive number eps and using the closed interval eps..1-eps. I don't know off the top of my head whether the limit of the solutions thus obtained will necessarily converge to the true solution as eps approaches 0; my guess is yes, but that there's no pratical way to compute how small eps needs to be.

Your procedures F and L are correct, although L's performance can be improved significantly by including option remember.

The procedure S is very wrong, and it's not entirely clear how to correct it. I think that you may mean this:

S:= proc(n::nonnegint)
local i;  
     for i from 1 by 1 to n do
          printf("%a\n", L(i)^(2)-5*F(i)^(2))
     end do  
end proc;

  1. There's no need to refer to S inside S.
  2. You were missing end do.
  3. Formatted printing is done with printf, not print.
  4. When using printf, you need to handle all spacing and line breaks. The \n is a "new line".

Your L can be shortened to

L:= proc(n::nonnegint)
option remember;
     if n = 0 then 2 elif n = 1 then 1 else L(n-1)+L(n-2) end if
end proc;

or

L:= proc(n::nonnegint)
option remember;
     if n = 0 then 2 elif n = 1 then 1 else thisproc(n-1)+thisproc(n-2) end if
end proc;

When you refer to a procedure from within it's own body, you can refer to it as thisproc. This is a tiny bit faster than referring to it by name, and it also means that you won't have to change the procedure's body at all if you decide to change the procedure's name.

 

Example:

P1:= plot3d(x*y, x= -1..1, y= -1..1):
P2:= plots:-pointplot3d([[0,0,0]], symbol= solidsphere, symbolsize= 32, color= red):
plots:-display([P1,P2]);

Why is your data organized in this strange way? Something tells me that you're making it more complicated than it needs to be. If each letter is to be paired with a number, then store them as pairs:

map(L-> {seq(L[2*k-1]=L[2*k], k= 1..nops(L)/2)}, data);

You see that when you use this more natural format, the sets are automatically sorted by letter. Unlike the other two Answers, mine will work on sublists of any length.

Your manual solution is faulty because the Laplace parameter s still appears in the equation. Since no gamma appears in your solution, I will assume that the coefficient gamma in the second original equation was supposed to be lambda.

These can be solved easily with Maple, like this:

restart:
Ps:= [P[0],P[1]]:
eqs:= Ps(s) =~ [1/(s+lambda)+mu*P[1](s)/(s+lambda), lambda*P[0](s)/(s+mu)]:
Ls:= solve(eqs, Ps(s))[];

Ps(t) =~ inttrans[invlaplace]~(rhs~(Ls), s, t);

If the step size is (b-a)/n, then including a and b, there are n+1, not n, evenly spaced points in the interval.

They can be generated like this (I take the name linspace from the equivalent Matlab function):

linspace:= (a,b,n)-> [seq(a+(b-a)*k/(n-1), k= 0..n-1)]:
linspace(0,1,50);

The procedure f can be applied to them like this:

f~(linspace(0,1,50));

While time is useful for testing a block of code (which should be in a single execution group), it is better to use CodeTools:-Usage to test a single procedure call. It returns more-detailed information. Example:

P:= `*`('randpoly(x, degree= 2^9, dense)' $ 2^7):
Prod:= CodeTools:-Usage(degree(expand(P)));
memory used=0.58GiB, alloc change=257.00MiB, cpu time=16.59s, real time=16.62s, gc time=0ns

     Prod:= 65535

By default, Usage returns directly the result of the computation and returns its usage report as a side effect. This default can be changed with Usage's output option.

This is a common source of bugs in Maple code: The code is expecting an expression sequence. If the sequence has only one member, then indexing the sequence becomes indexing of the member, which is not what was intended. The solution is to turn the expression sequence into a list by enclosing it in square brackets:

Answers:= [solve(equation that needs solving)];

This will work even if solve returns nothing. Now if Answers has only one solution, Answers[1] will return that solution in its entirety.

Terminology: Neither a list nor a sequence should be called a matrix, even when it can be double-indexed. And an expression sequence shouldn't be called a list, as in the title to this Question.

The number of elements of a list is returned by nops. This'll tell you what are the safe indices to use.

You should not use indexing to extract the a and b values from an individual solution. Use eval instead, as in

eval(a, Answers[1]), eval(b, Answers[1]);

Maple has tools to make the elimination of duplicates trivial and the elimination of near duplicates very easy.

Your code in both this Question and your previous Question seems to suggest that you have the mistaken belief that a data structure's indices need to be integers. This isn't true for a Maple table, which is what your x is. Its indices can be anything at all; they don't even need to be numbers. In the code below, I use the roots themselves as the indices, which automatically removes exact duplicates. Your attempt to use the integer indices i and j makes your code more complex.



restart:

Digits:= 15:
f:= z-> HankelH1(2,z):
x||(min,max,step):= (-50, 0, 0.2):
y||(min,max,step):= (-10, 0, 0.2):
eps:= 1e-5: #Root closeness tolerance

R:= table(): #To hold the roots

I only looked for roots in the third quadrant of the original rectangle that you used. This takes about an hour in Maple 16. It should be substantially quicker in Maple 18 and later because Hankel functions were added to evalhf.

st:= time():
for x from xmin by xstep to xmax do
     for y from ymin by ystep to ymax do
          r:= fsolve(f(z), z= x+I*y, z= xmin+I*ymin..xmax+I*ymax, complex);
          #If fsolve produces a numeric result, save it; otherwise discard.
          #Using the roots as the table's indices eliminates duplicates
          #automatically.
          if eval(r,1)::complexcons then
               R[r]:= [][] #Right side can be anything.
          end if
     end do
end do;
time()-st;

3658.161

(1)

Convert the table's indices to a list.

R1:= [indices(R, nolist)]:

nops(R1);

1629

(2)

Filter out any roots outside the desired rectangle.

R2:= remove(
     r-> xmin > Re(r) or Re(r) > xmax or ymin > Im(r) or Im(r) > ymax,
     R1
):

nops(R2);

20

(3)

The vast majority of the original roots were outside the rectangle!

 

Assume that roots closer than eps are actually the same root. Thus (r1,r2)-> abs(r1-r2) < eps is an equivalence relation. Select one member of each equivalence class.

R3:= map2(op, 1, [ListTools:-Categorize((r1,r2)-> abs(r1-r2) <= eps, R2)]):

nops(R3);

16

(4)

Check residuals.

max((abs@f)~(R3));

0.219000533285720e-13

(5)

R3;

[-14.7960226995499-.349557862119977*I, -49.4421659771279-.346839548441666*I, -17.9598588986496-.348595587873901*I, -1.31684116739175-.836104483291729*I, -30.5692124163275-.347269864797733*I, -11.6199893624061-.351427961366798*I, -40.0084502598130-.346979861585339*I, -8.41764502890611-.355893608952162*I, -46.2979989512446-.346876918993769*I, -27.4205845371005-.347439213958401*I, -21.1170212041082-.348034705428732*I, -43.1534565874607-.346922765296634*I, -36.8628610217597-.347052219103216*I, -24.2701281841205-.347679009437270*I, -33.7165254076117-.347145813344966*I, -5.13755909753331-.372212752744062*I]

(6)

The roots appear to lie along a "critical strip."

plot([Re,Im]~(R3), style= point, symbol= circle, symbolsize= 16);

 


plots:-complexplot3d(f(z), z= -50-.35*I..-20-.34*I);

 

plot(abs(f(X-.345*I)), X= -50..0, view= 0..0.5);

 

 



Download HankelRoots.mw

You have at least two errors with your final integral. One is a logic error, and the other is a syntax error. The logic error is that you use both x and y as variables and that you integrate with respect to x. I don't understand why you introduced x at all, and certainly the integration should be performed with respect to y.

The syntax error is that ay is being treated as a single variable. If you intend multiplication (as I'm sure you do), then you need to enter either a space between them (a y) or an explicit multiplication operator (a*y).

The issue with the decimal points is something called floating-point contagion: Once there's one decimal point in an expression, it tends to spread to other numbers in the expression automatically as operations are performed. In your case, the one decimal point comes from the 62.4. There are two ways you can prevent floating-point contagion. The first (and my preference) is to represent floating-point constants as symbols; the second is to represent them as exact rationals such as 624/10.

Use the seq command to create the lists sys and var, like this:

N:= 10:
sys:= [seq(galerkin_funcs[i], i= 1..N+4)]:
var:= [seq(w[i], i= 1..N)]:
(Kmat, Fmat):= LinearAlgebra:-GenerateMatrix(sys, var)

I did this without a ruler, just eyeballing it and restricting my coordinates to integers. You should be able to see how to make any numeric adjustments.

rectangle:= (P::[numeric,numeric], h::numeric, w::numeric)->
     plots:-polygonplot([P, P+~[w,0], P+~[w,h], P+~[0,h], P], args[4..])
:
     
plots:-display(
     [plot([[0,0], [10,0], [10,10], [0,10], [0,0]], thickness= 5, color= black),
      plots:-polygonplot([[0,10], [5,15], [10,10], [0,10]], color= red),
      rectangle([1,1], 3, 4, thickness= 2, color= white),
      rectangle([7,0], 4, 2, thickness= 2, color= blue),
      rectangle([1,7], 2, 2, thickness= 2, color= pink),
      rectangle([7,7], 2, 2, thickness= 2, color= pink),
      plots:-polygonplot([[8,12], [8,13], [9,13], [9,11], [8,12]], thickness= 2, color= grey)
     ], scaling= constrained, axes= none    

);

First 243 244 245 246 247 248 249 Last Page 245 of 395