Carl Love

Carl Love

28115 Reputation

25 Badges

13 years, 138 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

[Edit: It turns out that the first two of the three solution techniques presented in this Answer are irrelevant because the true cause of the problem was the OP's high setting of Digits (which was not revealed until much later in the discussion below). I just leave them here so that the rest of the thread makes sense.]

Here are three completely different techniques to solve this.

 Procedural form of input to Int and plot

Use the procedural form of input to Int and plot, which does not use an independent variable.

Also, your way of extracting the procedure from the dsolve solution is convoluted: There is no need to mention the independent variable. Ordinarily, what you did with the independent variables would make sense to Maple, although it is not necessary to do it. However, the procedures returned by dsolve(..., numeric, output= listprocedure) treat a symbolic independent variable in a slightly strange way: Instead of returning y(x), they return `y(x)`(x). This would be fine if the numeric procedure was named `y(x)`, but you've changed the name to u.

So, to correct the problem, do this:

# No change needed to dsolve.
So:= dsolve(
     [-1.2*diff(y(x), x$2)+0.8*y(x) = 2, y(0)=1, y(1)=0], y(x),
     'numeric', 'output' = listprocedure
):

# Simple extraction; no unapply needed.
u:= rhs(So[2]):

# No independent variables on the procedures or the range.
plot([u, u^2], 0..1);

# No independent variables, and note the capital I
evalf(Int(u, 0..1));

evalf(Int(u^2, 0..1));

 Or... Name the numeric procedure appropriately

You could just use `y(x)` instead of u as the target for the unapply, then make the next line u:= `y(x)`, and keep the rest of your code exactly as it was before:

`y(x)`:= unapply(rhs(So[2])(t), t):
u:= `y(x)`:

 Or... Make the integrals part of dsolve's output

You could simply include the appropriate functions in the original system passed to dsolve. It has often been argued here at MaplePrimes that this is the most numerically robust way to do it, due to dsolve's sophisticated error control. (However, I suspect that it is more computationally intensive because it requires the computation of the integrals at all the intermediary points.) This way involves a little bit of thought to rearrange the integrals into derivative form.

So:= dsolve([
     -1.2*diff(y(x), x$2)+0.8*y(x) = 2, y(0)=1, y(1)=0,

     diff(u(x),x) = y(x), u(0)=0,
     # So u(x) = int(y(t), t= 0..x).

     diff(u2(x),x) = y(x)^2, u2(0)=0
     # So u2(x) = int(y(t)^2, t= 0..x).

     ], [y(x), u(x), u2(x)],
     'numeric', 'output' = listprocedure
):

So(1);


The required integrals are the last two entries in the above list.

 

Do you mean that you want the display of the function to be J_0(k), i.e., a more standard notation than BesselJ? Then set

interface(typesetting= extended);

Then,

inttrans[fourier](piecewise(abs(x)<1,1/sqrt(1-x^2),0),x,k);

The typesetting level can also be changed from the Tools -> Options -> Display menu.

This is a bug. I'm amazed that it hasn't surfaced before now. Looking at the simple code for IntegrationTools:-Combine shows both the source of the bug and the workaround.

showstat(IntegrationTools:-Combine);

IntegrationTools:-Combine := proc(v)
local ints, x, X, input, output;
   1   ints := indets(v,'And(DefiniteIntegral,Not(MultipleIntegral))');
   2   x := IntegrationTools:-GetVariable(ints[1]);
   3   input := subsindets[':-flat'](v,'And(DefiniteIntegral,Not(MultipleIntegral))',t -> Change(t,GetVariable(t) = X));
   4   output := combine(input,'int');
   5   return subs(X = x,Change(output,X = x))
end proc

So, it is clear from lines 1 and 2 that the command can only work on definite integrals! For indefinite integrals, just use combine(..., int).

combine(Int(sin(x),x)+Int(cos(x),x), int);

plots:-polarplot([r, Pi/6, r=0..2]);

The upper limit for r is arbitrary; I chose 2 for no particular reason. Note that Pi is with a capital P.

This technique is akin to how one would plot a vertical line in rectangular (Cartesian) coordinates:

plot([3, y, y= -2..2]);

I'm not sure if command polarplot is available in Maple 13. If it isn't, I'll come up with something else for you.

Perhaps they meant semicolons, which are used to separate rows:

< 6,8 ; -3,9 ; 4,14 >;

But that's row-major order.

By the way, using the angle-bracket operators uses significantly less processor time than using Matrix or Vector, so I'd recommend using them (or the built-in rtable) in programs whenever possible. Often, it is useful or necessary to use the prefix forms `<,>` and `<|>`. (There is no `<;>`.) For example, to convert a row-major listlist LL into a Matrix, you could do

<`<|>`~(LL)[]>;

But that example is just to illustrate the power of the prefix notation, because instead I'd probably use

rtable(LL, subtype= Matrix);

If your answer to the first question in my Reply above is yes, then this routine will do it:

SeparateLastTerms:= E->
     evalindets(
          E, specfunc(anything,{Sum,sum}),
          S-> eval(op(1,S), op([2,1],S)= op([2,2,2],S)) + subsop([2,2,2]= op([2,2,2],S) - 1, S)
     )
;

Example of use:

E:= Sum(A[k], k= 2..n+1) + sum(f(j), j= 3..m):
SeparateLastTerms(E);

Change the line

evalm(Matrix(n,n,symbol=C))

to

Matrix(n,n, (i,j)-> C[i,j])

The symbol=C option does not work when C has already been assigned values.

There is no need---ever---to use evalm with capital-M Matrices, which are the only type of Matrices that you should be using. Likewise, there is no need for with(linalg).

You wrote:

I understand in Maple one uses the back quote key (or rather the apostrophe, 0X27) to prevent one time evaluation of expression.

The apostrophe is not the same as the back quote! The aposthrope is ASCII character 39 (or 0X27 as you say); the back quote is `, ASCII character 96 (or 0X60). It is important in Maple not to mix these up! They each have very specific uses. The aposthrophe delays evaluation, and the back quote turns arbitrary groups of characters into names.

This might indicate that the front end parser did this simplification before the main evaluator got hold of it, so it was too late?

Yes, that is essentially correct. We refer to it as "automatic simplification." Yes, these cannot be delayed with quotes.

Either way, how would one make Maple return 16/4 when the input is '16/4'?

In Maple 18, use the InertForm package:

InertForm:-Display(InertForm:-Parse("16/4"));

or

InertForm:-Display(`%/`(16,4));

(Note the use of back quotes in the last command.)

You need to choose a branch of the absolute value. Each branch gives a different solution. When there are multiple solutions possible for the highest derivative, numeric dsolve will not choose one of them for you. So, solve the systems y*y''' = -1 and y*y''' = 1 separately and look at the different plots.

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]]);

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