Carl Love

Carl Love

28025 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

There are two steps to convert a list of points into piecewise form. First, apply the procedure JoggedLine (which can be found in this post by Kitonum, the procedure's author) to the list of points. This will give you an expression in absolute value form. Then apply convert(..., piecewise) to that expression. Example:

JoggedLine([[0,2],[3,-2],[7,6]]);

convert(%, piecewise);

 

Adri's Answer suffers from a common flaw in Maple code: It builds a sequence or list by iteratively appending each new element onto the end of an existing sequence or list. This is very inefficient in Maple because on each iteration it creates a new sequence or list. Thus the time complexity of this operation is O(n^2) where n is the number of items in the sequence or list.

There are two main alternatives: (1) Build with a seq command, or (2) build using a table or dynamic Vector or Array and convert to a sequence or list as the last step. The seq option is often not possible, but it is in this case. The are also several other ways to simplify the code.


restart:

S:= proc(a, b, s)

local x, y;
     Array([seq([y, min(select(`>`, [fsolve(f(x,y))], 0))], y= a..min(b,1), s)])

end proc:

 

Example of use:

f:= unapply(randpoly([x,y]), x,y);

proc (x, y) options operator, arrow; -92*x^3*y^2+6*x^2*y^3+74*x*y^4+23*x^4+75*x*y^3-50*x end proc

R:= S(.2, 2, .1);

R := Matrix(9, 2, {(1, 1) = .2, (1, 2) = 1.34423469254842, (2, 1) = .3, (2, 2) = 1.40236242650571, (3, 1) = .4, (3, 2) = 1.48386699558845, (4, 1) = .5, (4, 2) = 1.59489085298900, (5, 1) = .6, (5, 2) = 1.75109495862678, (6, 1) = .7, (6, 2) = 1.98664891453434, (7, 1) = .8, (7, 2) = .697772717030054, (8, 1) = .9, (8, 2) = 1.08004688099483, (9, 1) = 1.0, (9, 2) = 1.31587196629247})

 

``


Download Array_build.mw


restart:

The procedure:

Solve:= (A::Matrix(square), B::Matrix(square), C::Matrix(square))->
     LinearAlgebra:-LinearSolve(A+B, C):

Example of use:

(A,B,C):= ('LinearAlgebra:-RandomMatrix(2,2)' $ 3);

A, B, C := Matrix(2, 2, {(1, 1) = 44, (1, 2) = -31, (2, 1) = 92, (2, 2) = 67}), Matrix(2, 2, {(1, 1) = 8, (1, 2) = 99, (2, 1) = 69, (2, 2) = 29}), Matrix(2, 2, {(1, 1) = -32, (1, 2) = -4, (2, 1) = -74, (2, 2) = 27})

x:= Solve(A,B,C);

x := Matrix(2, 2, {(1, 1) = -490/1489, (1, 2) = 555/1489, (2, 1) = -326/1489, (2, 2) = -512/1489})

Verification of correctness:

A.x + B.x;

Matrix(2, 2, {(1, 1) = -32, (1, 2) = -4, (2, 1) = -74, (2, 2) = 27})

LinearAlgebra:-Equal(%, C);

true

``


Download Matrix_solve.mw

You can trap the error with the try statement. Here's an example. It is not necessary to use parse; any code can come after the try. The parse statement is useful if you want to set up your 50 problems as strings of code to be evaluated.

restart:
Problems:= Vector(50, k-> sprintf("50/(25-%d)", k)):
Answers:= Vector(50):
for k to 50 do
     try
          Answers[k]:= parse(Problems[k])
     catch:
          Answers[k]:= sprintf("k = %d: %s", k, lasterror);
          print('k'= k, "Error found")
     end try
end do;

Answers[25];

"k = 25: numeric exception: division by zero"

Answers[50];

−2

See ?try .

For your second question T(..., ...) is called "programmer indexing". See ?rtable_indexing .

The following is a hack that temporarily patches the deficiency/bug in PDEtools:-dchange that Alejandro pointed out. It gets you what you aked for: the application of dchange, but with unexpanded/unevaluated derivatives.

 

restart:

Reynolds:= Diff(p(x)*h(x)^3/(12*mu)
     *Diff(p(x),x),x)-u(x)/2*Diff(p(x)*h(x),x)
     +Diff(p(x)*h(x)^3/(12*mu)*Diff(p(x),z),z)
     =Diff(p(x)*h(x),t);

Diff((1/12)*p(x)*h(x)^3*(Diff(p(x), x))/mu, x)-(1/2)*u(x)*(Diff(p(x)*h(x), x))+Diff((1/12)*p(x)*h(x)^3*(Diff(p(x), z))/mu, z) = Diff(p(x)*h(x), t)

Start of hack.

unprotect(diff);  sav_diff:= eval(diff):  diff:= fiff: #Nonsense name.

ans1:= PDEtools:-dchange({p(x)=P(X)*Pa,x=Lx*X,h(x)=H(X)*h2},Reynolds,{P,h,X,u,H});

fiff(fiff((1/12)*P(X)*Pa*H(X)^3*h2^3*fiff(fiff(P(X)*Pa, X)/fiff(Lx*X, X), [])/mu, X)/fiff(Lx*X, X), [])-(1/2)*u(X, Lx)*fiff(fiff(P(X)*Pa*H(X)*h2, X)/fiff(Lx*X, X), [])+fiff((1/12)*P(X)*Pa*H(X)^3*h2^3*fiff(P(X)*Pa, [z])/mu, [z]) = fiff(P(X)*Pa*H(X)*h2, [t])

diff:= eval(sav_diff):  protect(diff):

ans2:= eval(ans1, [fiff= ((f,x)-> `if`(x=[], f, Diff(f,x)))]);

(Diff((1/12)*P(X)*Pa*H(X)^3*h2^3*(Diff(P(X)*Pa, X))/(mu*(Diff(Lx*X, X))), X))/(Diff(Lx*X, X))-(1/2)*u(X, Lx)*(Diff(P(X)*Pa*H(X)*h2, X))/(Diff(Lx*X, X))+Diff((1/12)*P(X)*Pa*H(X)^3*h2^3*(Diff(P(X)*Pa, [z]))/mu, [z]) = Diff(P(X)*Pa*H(X)*h2, [t])

fiff:= eval(diff):

End of hack.

ans2 is the desired answer with unevaluated derivatives. The rest of this worksheet is simply a proof of the correctness of ans2.

save Reynolds, ans2, "C:/Users/Carl/desktop/dchange_hack.mpl";


restart:

read "C:/Users/Carl/desktop/dchange_hack.mpl":

Use regular PDEtools:-dchange.

ans3:= PDEtools:-dchange({p(x)=P(X)*Pa,x=Lx*X,h(x)=H(X)*h2},Reynolds,{P,h,X,u,H});

((1/12)*(diff(P(X), X))^2*Pa^2*H(X)^3*h2^3/(mu*Lx)+(1/4)*P(X)*Pa^2*H(X)^2*h2^3*(diff(P(X), X))*(diff(H(X), X))/(mu*Lx)+(1/12)*P(X)*Pa^2*H(X)^3*h2^3*(diff(diff(P(X), X), X))/(mu*Lx))/Lx-(1/2)*u(X, Lx)*((diff(P(X), X))*Pa*H(X)*h2+P(X)*Pa*(diff(H(X), X))*h2)/Lx = 0

value(ans2);

((1/12)*(diff(P(X), X))^2*Pa^2*H(X)^3*h2^3/(mu*Lx)+(1/4)*P(X)*Pa^2*H(X)^2*h2^3*(diff(P(X), X))*(diff(H(X), X))/(mu*Lx)+(1/12)*P(X)*Pa^2*H(X)^3*h2^3*(diff(diff(P(X), X), X))/(mu*Lx))/Lx-(1/2)*u(X, Lx)*((diff(P(X), X))*Pa*H(X)*h2+P(X)*Pa*(diff(H(X), X))*h2)/Lx = 0

simplify(lhs(ans3) - lhs(%));

0

 

``

 

Download dchange_hack.mw

I have a feeling that there's a more efficient algorithm to generate a round-robin tournament schedule than what I have below. I end up throwing out 2/3 of the partitions. Hopefully a more-skilled combinatorist can comment on the efficiency. Nonetheless, it does get you the answer that you need.

Other interesting questions are:

  1. How do you count how many essentially distinct round-robin tournament schedules there are?
  2. Is it possible to generate them efficiently with an iterator?
  3. How do the answers change when the number of teams is odd? (The code below will not work when the number of teams is odd.)

In the code below, remember to use compile= false as an option to SetPartitions if you don't have a Maple compiler.

restart:
T:= {}:
night:= 0:
for P in Iterator:-SetPartitions({A,B,C,D,E,F}, [[2,3]]) do
     p:= {P[]}:
     if p intersect T = {} then
          T:= T union p;
          night:= night+1:
          printf("\nNight %d:", night);
          for s in p do
               printf(" %a vs. %a ", s[])
          end do
      end if
end do:

Night 1: A vs. B  C vs. D  E vs. F
Night 2: A vs. C  B vs. E  D vs. F
Night 3: A vs. D  B vs. F  C vs. E
Night 4: A vs. F  B vs. C  D vs. E
Night 5: A vs. E  B vs. D  C vs. F

Solve symbolically first. Use unapply to turn the symbolic solution into a procedure with paramter a. Then map that procedure over Vector a. Doing it this way, I accomplished the entire job in 0.2 seconds.

restart:
n:= 10^4:
eq1:={b = a^2, c = b^3/2, d = c^(1/2)*4 + b^2}: #No quote marks on a!
st:= time():
Sol:= unapply(solve(eq1, {b,c,d}), a);

a:= Vector(n, i-> i):
ans:= map(Sol, a):
time()-st;
                            
0.188


restart:

u:= -6*a+b+c:  v:= 3*a-4*b+c:

expand(u*v)/sqrt(expand(u*u)*expand(v*v));

(-18*a^2+27*a*b-3*a*c-4*b^2-3*b*c+c^2)/((36*a^2-12*a*b-12*a*c+b^2+2*b*c+c^2)*(9*a^2-24*a*b+6*a*c+16*b^2-8*b*c+c^2))^(1/2)

eval(%, [a*b=1, a*c=1, b*c=1, a^2=1, b^2=16, c^2=16]):

arccos(%);

Pi-arccos((3/782)*11730^(1/2))

evalf(%);

1.99928084710157

evalf(convert(%, degrees));

114.550354600260*degrees

 

``


Download AngleBetween.mw

Let n represent the number of degrees of freedom. Then

X:= Statistics:-RandomVariable(StudentT(n)):
T:= unapply(Statistics:-CDF(X,x), n,x);

plot(T(1,x), x= -3..3);

evalf(T(1,1));
                       0.750000000000000

Please ask any further questions about logic problems in a new thread.

Here is how to obtain a complete solution to the logic problem, and to prove its uniqueness, using my LogicProblem package (available at the Maple Applications Center).

Solved 2014-Feb-10 from Puerto Viejo de Talamanca, Limon, Costa Rica.

restart:

Statement of the logic problem:

Who Works Where?

Alex, Betty, Carol, Dan, Earl, Fay, George and Harry are eight employees of an organization

 

Constraint 0: They work in three departments: Personnel, Administration and Marketing with not more than three of them in any department.

 

Each of them has a different choice of sports from Football, Cricket, Volleyball, Badminton, Lawn Tennis, Basketball, Hockey and Table Tennis not necessarily in the same order.

 

Constraints:

1. Dan works in Administration and does not like either Football or Cricket.

2. Fay works in Personnel with only Alex who likes Table Tennis.

3. Earl and Harry do not work in the same department as Dan.

4. Carol likes Hockey and does not work in Marketing.

5. George does not work in Administration and does not like either Cricket or Badminton.

6. One of those who work in Administration likes Football.

7. The one who likes Volleyball works in Personnel.

8. None of those who work in Administration likes either Badminton or Lawn Tennis.

9. Harry does not like Cricket.

 

Questions:

1. Who are the employees who work in the Administration Department?

2. In which Department does Earl work?

 

Vars:= [Name, Dept, Sport]:

Name:= [alex, betty, carol, dan, earl, fay, george, harry]:

Constraints 0 & 2 together imply that the departments are divided thus:

person:= {person||(1..2)}:  admin:= {admin||(1..3)}:  market:= {market||(1..3)}:

Dept:= [person[], admin[], market[]]:

Sport:= [football, cricket, volleyball, badminton,
      lawn_tennis, basketball, hockey, table_tennis]:

Con1:= dan = admin1, (dan <>~ {football, cricket})[]:

Con2:= fay = person1, alex = person2, alex = table_tennis:

Con3:= (earl <>~ admin)[], (harry <>~ admin)[]:

Con4:= carol = hockey, (carol <>~ market)[]:

Con5:= (george <>~ admin union {cricket, badminton})[]:

Con6:= (football <>~ market union person)[]:

Con7:= (volleyball <>~ market union admin)[]:
Con8:= (badminton <>~ admin)[], (lawn_tennis <>~ admin)[]:

Con9:= harry <> cricket:

read "C:/Users/Carl/desktop/logic_problems.mpl":

Work:= LogicProblem(Vars):

with(Work);

Warning, Work is not a correctly formed package - option `package' is missing

[`&!!`, `&-`, `&<`, `&>`, `&?`, `&G`, `&Soln`, AutoGuess, CPV, CollectStats, ConstNum, Consts, ConstsInV, DifferentBlock, Equation, FreeGuess, GoBack, Guess, InternalRep, IsComplete, IsUnique, NC, NV, PrintConst, Quiet, Reinitialize, SameBlock, Satisfy, Separated, UniquenessProof, VarNum, VarNumC, X_O]

Satisfy([Con||(1..9)]);

NULL

`Incomplete solution; need more constraints to complete:`

Matrix(8, 3, {(1, 1) = alex, (1, 2) = person2, (1, 3) = table_tennis, (2, 1) = betty, (2, 2) = _, (2, 3) = _, (3, 1) = carol, (3, 2) = _, (3, 3) = hockey, (4, 1) = dan, (4, 2) = admin1, (4, 3) = basketball, (5, 1) = earl, (5, 2) = _, (5, 3) = _, (6, 1) = fay, (6, 2) = person1, (6, 3) = volleyball, (7, 1) = george, (7, 2) = _, (7, 3) = _, (8, 1) = harry, (8, 2) = _, (8, 3) = _})

But the solution is complete enough to answer the two questions posed: Who are the employees that work in the Administration Dept.? In which department does Earl work?

Name &? Dept;

Matrix(9, 9, {(1, 1) = ` `, (1, 2) = alex, (1, 3) = betty, (1, 4) = carol, (1, 5) = dan, (1, 6) = earl, (1, 7) = fay, (1, 8) = george, (1, 9) = harry, (2, 1) = person1, (2, 2) = X, (2, 3) = X, (2, 4) = X, (2, 5) = X, (2, 6) = X, (2, 7) = O, (2, 8) = X, (2, 9) = X, (3, 1) = person2, (3, 2) = O, (3, 3) = X, (3, 4) = X, (3, 5) = X, (3, 6) = X, (3, 7) = X, (3, 8) = X, (3, 9) = X, (4, 1) = admin1, (4, 2) = X, (4, 3) = X, (4, 4) = X, (4, 5) = O, (4, 6) = X, (4, 7) = X, (4, 8) = X, (4, 9) = X, (5, 1) = admin2, (5, 2) = X, (5, 3) = _, (5, 4) = _, (5, 5) = X, (5, 6) = X, (5, 7) = X, (5, 8) = X, (5, 9) = X, (6, 1) = admin3, (6, 2) = X, (6, 3) = _, (6, 4) = _, (6, 5) = X, (6, 6) = X, (6, 7) = X, (6, 8) = X, (6, 9) = X, (7, 1) = market1, (7, 2) = X, (7, 3) = _, (7, 4) = X, (7, 5) = X, (7, 6) = _, (7, 7) = X, (7, 8) = _, (7, 9) = _, (8, 1) = market2, (8, 2) = X, (8, 3) = _, (8, 4) = X, (8, 5) = X, (8, 6) = _, (8, 7) = X, (8, 8) = _, (8, 9) = _, (9, 1) = market3, (9, 2) = X, (9, 3) = _, (9, 4) = X, (9, 5) = X, (9, 6) = _, (9, 7) = X, (9, 8) = _, (9, 9) = _})

The Xs in the table represent excluded possibilities, the Os represent known equivalences, and the _s are the equivalences that have been neither proven nor disproven.

From this table we see that Dan, Betty, and Carol work in Administration and that Earl works in Marketing.

 

It is possible to get to a unique solution with the given constraints. Add the information that we just found out.

Satisfy([betty = admin2, earl = market1, george = market2]);

NULL

`Unique solution:`

Matrix(8, 3, {(1, 1) = alex, (1, 2) = person2, (1, 3) = table_tennis, (2, 1) = betty, (2, 2) = admin2, (2, 3) = football, (3, 1) = carol, (3, 2) = admin3, (3, 3) = hockey, (4, 1) = dan, (4, 2) = admin1, (4, 3) = basketball, (5, 1) = earl, (5, 2) = market1, (5, 3) = cricket, (6, 1) = fay, (6, 2) = person1, (6, 3) = volleyball, (7, 1) = george, (7, 2) = market2, (7, 3) = lawn_tennis, (8, 1) = harry, (8, 2) = market3, (8, 3) = badminton})

 

 

Download WhoWorksWhere.mw

You need to put the if statement inside the procedure, like this:

MARK:= (a,b,c,d)->
    if a.d-b.c = 0 then
         [0,0,0,0]
    else
         [d/(a.d-b.c), -b/(a.d-b.c), -c/(a.d-b.c), a/(a.d-b.c)]
    fi;

The print statement is not needed. Indeed, print is never an appropriate way for a procedure to return its output.

restart:
`mod/ReduceExponents`:= (p::polynom, n::posint)->
     evalindets(p, '`^`'(name, posint), T-> op(1,T)^(op(2,T) mod n)):
p:= x-> 1+x^3+x^17+x^19:
ReduceExponents(p(x)) mod 17;

ArrayInterpolation will not accept exact integer values; they must be converted to floating-point values with evalf. Also, do not use with inside a procedure; use uses instead. And I made your variables local. A will be the return value of the procedure.

InterProc := proc ()
uses CurveFitting;
local i, A, variable1, variable2;
     A := Matrix(10, 4);
     variable1 := evalf([0, 100]);
     variable2 := evalf([12, 20]);

     for i from 1 to 10 do
          A(i, 1) := evalf(3+i);
          A(i, 2) := 2*i+A(i, 1);
          A(i, 3) := A(i, 2)-1;
          A(i, 4) := ArrayInterpolation(variable1, variable2, A(i, 1));
     end do;
end proc;

Three errors I see:

  1. "if" must be with a lowercase "i".
  2. Need a semicolon at the end of the first line.
  3. The denominators need to be in parentheses: (a*d - b*c).

You should make your code into a procedure, like below. There is no need for the print command.

 

restart:

Inverse:= proc(M::Matrix(2,2))
local a,b,c,d,D;
    (a,b,c,d):= convert(M,list)[];
    D:= a*d-b*c;
    if D=0 then
         error "Matrix is not invertible."
    else
         < < d, -b > | < -c, a > > / D
    end if
end proc:

M:= < < a, b > | < c , d > >;

M := Matrix(2, 2, {(1, 1) = a, (1, 2) = c, (2, 1) = b, (2, 2) = d})

Inverse(M);

Matrix(2, 2, {(1, 1) = d/(a*d-b*c), (1, 2) = -c/(a*d-b*c), (2, 1) = -b/(a*d-b*c), (2, 2) = a/(a*d-b*c)})

M.Inverse(M);

Matrix(2, 2, {(1, 1) = a*d/(a*d-b*c)-b*c/(a*d-b*c), (1, 2) = 0, (2, 1) = 0, (2, 2) = a*d/(a*d-b*c)-b*c/(a*d-b*c)})

simplify(%);

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})

M:= Matrix(2,2);

M := Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0})

Inverse(M);

Error, (in Inverse) Matrix is not invertible.

M:= LinearAlgebra:-RandomMatrix(2,2);

M := Matrix(2, 2, {(1, 1) = 44, (1, 2) = -31, (2, 1) = 92, (2, 2) = 67})

Inverse(M);

Matrix(2, 2, {(1, 1) = 67/5800, (1, 2) = 31/5800, (2, 1) = -23/1450, (2, 2) = 11/1450})

M.Inverse(M);

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = 1})

 

 

Download Inverse.mw

 

Use the function form on input, delta(x) and Delta(x), and define the following procedures, perhaps in your initialization file:

`diff/delta`:= x-> delta(x):
`diff/Delta`:= x-> Delta(x):
`print/delta`:= x-> delta*x:
`print/Delta`:= x-> Delta*x:

Now the function form will display as the multiplication form, it won't be affected by taking the derivative, and simplify will only see the function form and won't change it.

First 320 321 322 323 324 325 326 Last Page 322 of 395