Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The command is ListTools:-Rotate.

Your code may be syntactically correct in the sense that the parser accepts it, but I would not call it syntactically correct in the broader sense. Your call to DirectSearch:-Search contains the constraint

dist(a, b, x_1, y_1, z_1, x_2, y_2, z_2, phi_1, psi_1, theta_1, phi_2, psi_2, theta_2)[1] >= 0.1e-1

This forces a call to dist with symbolic arguments, which, of course, is nearly impossible to compute, although DirectSearch:-GlobalOptima apparently tries. Note that this call to DirectSearch:-GlobalOptima occurs before the call to DirectSearch:-Search because it happens when the arguments are being evaluated before being passed. To correct the problem, change the constraint to procedural form:

(a, b, x_1, y_1, z_1, x_2, y_2, z_2, phi_1, psi_1, theta_1, phi_2, psi_2, theta_2)->
     evalb(
          dist(a, b, x_1, y_1, z_1, x_2, y_2, z_2, phi_1, psi_1, theta_1, phi_2, psi_2, theta_2)[1]
          >= 0.01
     )

Then dist will only be called with numeric arguments.

 

The derivative of almost any numeric or symbolic procedure can be done with D; it works on a much wider range of procedures than does codegen[GRADIENT].

ODE:= diff(y(x),x$2) - 1.0325*diff(y(x),x) + 1.36*y(x) = sin(2*x):
S1:= dsolve([ODE, y(0)=0, y(1)=1], numeric, y(x), 'output'= listprocedure):
D2:= D(eval(diff(y(x),x), S1)):

plot(D2(t), t= 0..1, thickness= 3);

However, I suspect (but I am not sure) that the second derivative constructed from the original ODE is more numerically stable than the one produced by D. That is because in this case does not differentiate the code of the procedure (it cannot understand the complex code output by dsolve), but rather it constructs a template for numeric differentiation. In the case of simpler procedures, D actually differentiates the procedures' code just like codegen[GRADIENT] does.

 

It can be easily done with map2 and the elementwise operator ~; no explicit looping is required.

The package linalg is old and deprecated. Replace calls to its commands with calls to the equivalent commands from the packages VectorCalculus or LinearAlgebra.

In my code below, note the correct way to construct the Jacobian procedure Df with unapply. The way that you were using, with ->, was recomputing the Jacobian symbolically for every pair of numeric arguments---very inefficient.

restart:

kappa:= 2:  g:= 1:

sys:= [u*(1-u/kappa)-u*v, g*(u-1)*v]:

Vars:= [u,v]:

#Compute steady states:
SSs:= map2(eval, Vars, [solve(sys)]);

[[0, 0], [2, 0], [1, 1/2]]

Df:= unapply(VectorCalculus:-Jacobian(sys, Vars), Vars):

Df(Vars[]);

Matrix(2, 2, {(1, 1) = 1-u-v, (1, 2) = -u, (2, 1) = v, (2, 2) = u-1})

#All the Jacobian matrices:
Mats:= (Df@op)~(SSs);

Mats := [Matrix(2, 2, {(1, 1) = 1, (1, 2) = 0, (2, 1) = 0, (2, 2) = -1}), Matrix(2, 2, {(1, 1) = -1, (1, 2) = -2, (2, 1) = 0, (2, 2) = 1}), Matrix(2, 2, {(1, 1) = -1/2, (1, 2) = -1, (2, 1) = 1/2, (2, 2) = 0})]

#All the eigenvalues:
LinearAlgebra:-Eigenvalues~(Mats);

[Vector(2, {(1) = 1, (2) = -1}), Vector(2, {(1) = 1, (2) = -1}), Vector(2, {(1) = -1/4+((1/4)*I)*sqrt(7), (2) = -1/4-((1/4)*I)*sqrt(7)})]

 

 

Download Predator_Prey.mw

 

Here's a somewhat automated way to get a numeric procedure for the second derivative using the solution returned by dsolve and the original ODE.

ODE:= diff(y(x),x$2) - 1.0325*diff(y(x),x) + 1.36*y(x) = sin(2*x):
S1:= dsolve([ODE, y(0)=0, y(1)=1], numeric, y(x), 'output'= listprocedure):
D2:= unapply(subsindets(eval(solve(ODE, diff(y(x),x$2)), S1[2..-1]), procedure, P-> 'P'(x)), x):

plot(D2(x), x= 0..1, thickness= 3);

Note that the x in the plot command is superfluous. I could just as well have used plot(D2, 0..1). I included the x to emulate your coding style from the original plot.

This is known as the one-dimensional bin-packing problem. See the Wikipedia article "Bin-packing problem".

Here's an even more programmatic solution. In particular, I have automated the count of the maximum number of files per disk. And I put the output in a more convenient format. The algorithm is the same as Kitonum's. This algorithm is probably not very good if very many of the smallest file can fit on a disk or if there are very many different file sizes.

BinPacking1D:= proc(
     Size::positive, #Size of bins or disks
     #List of pairs: [object size, number of objects of that size]:
     Objects::list([positive, posint])
)
uses C= combinat, LT= ListTools;
local
     Szs:= map2(op, 1, Objects), #Sizes
     Cts:= table((`=`@op) ~ (Objects)), #Counts
     All:= (`$`@op) ~ (Objects),
     m, size, g,
     S:= select(g-> `+`(g[]) <= Size, [seq(C:-choose(All, m)[], m= 1..trunc(Size/min(Szs)))]),
     V:= table([seq(g= nprintf("%a", g), g= S)]) #Decision variables    
;
     (e-> parse(lhs(e))=rhs(e)) ~
          (select(
               e-> rhs(e)<>0,
               Optimization:-LPSolve(
                    `+`(entries(V))[],
                    [seq(add(V[g]*LT:-Occurrences(size, g), g= S) = Cts[size], size= Szs)],
                    assume= [integer, nonnegative]
               )[2]
          ))
end proc:

BinPacking1D(1.44, [[0.8, 3], [0.7, 12], [0.4, 15]]);

You have three errors:

1. You cannot use square brackets or curly braces to substitute for parentheses. Only parentheses can be used for algebraic grouping. The other brackets and braces have their own meanings.

2. The letter e means nothing to Maple on input; you need to use the function exp. So, e^x becomes exp(x).

3. Remove the "=energy" part of the equation; use only the left side. The plot3d command plots expressions or procedures, not equations. 

The unknowns are _C1, ..., _C4, not C1, ..., C4. If you include the underscores, the solve command will quickly return your solutions.

solution:= solve({t11 = 0, t22 = 0, t3 = -B, t4 = -A}, {_C1, _C2, _C3, _C4});

In the 2D input, which is what you are using, I believe that you must enter the underscores by typing \_ (backslash underscore).

You have a bug in the line

M_finish:= (X_1,X_2,X_3,X_4,X_5,X_6,z)->
    
M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X_5+M_6*X_6+M_F:

I'll use a simplified example to explain the bug:

restart:
M:= 2*z:
f:= z-> M:
f(1);

f(2);

You see, the z that is inside M is the global z and not the z that is the parameter of f. The correct way to construct the procedure f is with unapply:

f:= unapply(M,z);

f(1);

      2


 

I can't tell what you are doing wrong, but I suspect that it has something to do with using 2d input, which I distrust. It seems to be ignoring some of your commands as if they were commented out.

Using 1d input, and doing essentially the same thing as you (I added variable initial conditions):

restart:
eq1:= diff(V(t),t) = alpha-beta:
Q:= t-> V(t)*C(t):
eq2:= diff(Q(t),t) = G - K*C(t):
Sol:= dsolve({eq1, eq2, V(0)= V0, C(0)= C0});


Maple will not return an answer for a numeric integral unless it has some confidence that the error is less than what is specified by Digits (by default, it is 10 significant digits). Set the Digits environment variable higher and set the digits option to Int lower.

restart:
Digits:= 30:
#I distributed the sqrt(x) in the inner integral.
Inner:= Int(
     0.579160e-1*Ei(1., 0.500000e-4*x)+(1072.23*(.999950-1.*exp(-0.500000e-4*x)))/x,
     x = 1. .. eta,
     digits= 15
):
Int(exp(18.1818*Inner-9.10000*eta)/eta, eta = 1. .. 100., digits= 5):
evalf(%);

     0.00046645

I see that you have six digits of precision in your constants. I didn't have the patience to ge the answer to six digits, but I think it will finish. I would take it out further than six digits to be confident of the digits.

I found a closed-form formula for phi[j](x,y).



restart:

phi[1]:= (x,y)-> exp(x-y):

phi[j]:= (x,y)-> int(exp(x-s)*phi[j-1](s,y), s= y..x);

proc (x, y) options operator, arrow; int(exp(x-s)*phi[j-1](s, y), s = y .. x) end proc

(1)

Phi:= (j,x,y)-> exp(x-y)*(x-y)^(j-1)/(j-1)!:

 

Proposition: For all positive integers j,

'phi[j](x,y)' = Phi(j,x,y);

phi[j](x, y) = exp(x-y)*(x-y)^(j-1)/factorial(j-1)

(2)

Proof: Use induction on j.

Base case: j=1:

simplify(Phi(1,x,y) - phi[1](x,y));

0

(3)

Inductive step:

subs(phi[j-1](s,y)= Phi(j-1,s,y), phi[j](x,y) - Phi(j,x,y));

int(exp(x-s)*exp(s-y)*(s-y)^(j-2)/factorial(j-2), s = y .. x)-exp(x-y)*(x-y)^(j-1)/factorial(j-1)

(4)

simplify(%) assuming j::posint;

0

(5)


Download Induction.mw

Here's another way: We make a single recursive procedure that reads its own index (whatever's in the square brackets after the procedure name). No loop is needed.

phi:= proc(x,y)
local s, j:= op(1,procname);
     `if`(j=1, exp(x-y), int(exp(x-s)*thisproc[jj-1](s,y), s= y..x))
end proc:


phi[1](x,y);

phi[2](x,y);

phi[3](x,y);

 

Use plots:-odeplot. For example, assuming your independent variable is t:

plots:-odeplot(dsn, [t, Y[1](t)], t= 0..10);

You have three syntax errors. I'm surprised that you did not get an error message.

The first two errors are that you used both proc(x) ... and x-> ... to define procedures. Every procedure must begin with one of these, but never both.

The third error is that you didn't declare the variable of integration. Given a procedure f, its integral can be expressed these two ways: int(f, a..b) or int(f(x), x= a..b). (There a few other syntaxes also.)

So here's your code, corrected. I pulled the procedure definitions out of the loop and simplified a complicated logical expression by using xor. And I removed needless repetition from the sum of integrals expression.

PMmmf_func:= x->
     `if`(x >= rotorshift xor trunc((x-rotorshift)/tau[p])::odd, 1, -1)*H[c]*l[m]:

B[g]:= x-> 1000*mu[0]*(PMmmf_func(x)+MMF_func(x))/d[eff, stator](x):

for i from 0 to 42 do
     rotorshift:= evalf[6](2*(1/180* i *Pi*(1/2*OO[gap]))/N[m]-1/2*(tau[p]-tau[s]));
     Flux(i+1):=
          add(Int(B[g](x), x= map(`*`, R, tau[s])), R= [0..1, 3..4, 6..7, 9..10])*1e-3*L[ro]    
end do:

This will return the integrals in unevaluated Int form. If you want numeric results, apply evalf to these.

 

First 281 282 283 284 285 286 287 Last Page 283 of 395