Carl Love

Carl Love

28095 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

Please don't check off all the "Products" categories in the message header. Your Question is about Maple itself, not about any other product.

Getting a result with a RootOf isn't an error. It's an indication that there's no simpler way of expressing the exact roots of your equation. If you want a decimal approximation of the real root, use fsolve(EQ). If you want decimal approximations of all nine complex roots, use fsolve(EQ, complex).

Using the view option will stabilize the viewing window, thus making it appear that the point moves.

Explore(plots:-pointplot(a*[1, 1], view= [-10..10, -10..10]), a= -10..10);

I don't use context menus, so I don't know if the above is accessible via one.

I've only analyzed the case phi[4]*S[3]; I guess that the analysis for the other cases is similar. Your integrand has a singularity at x=1. ApproximateInt doesn't handle improper integrals.

If you don't provide any extra options, then ApproximateInt uses a midpoint Riemann method with only 11 evenly spaced evaluation points. ApproximateInt is primarily intended for educational illustration purposes. There is no inherent error control. You would have to decide on the value of partition required to achieve a certain accuracy.

Regarding Int/int's epsilon option: If you set epsilon to .5*10^(1-d), then the result will have d digits of accuracy if Maple can complete the calculation to that accuracy without using too many extra guard digits. If Maple can't complete the calculation thus, the integral is returned unevaluated. In that case, increasing the value of Digits while keeping epsilon the same may help the computation to complete.

The k in your expression is used as a function symbol, not a coefficient. A number can be used as a function that is a constant function that returns the number. For example, 3(x) is simply 3. There is no d in your expression, so I don't know why you expect something to happen when you substitute for d.

To enter boundary or initial conditions that have derivatives, use D form rather than diff form.

Maple's numeric dsolve won't allow infinity as a boundary point. You'll just have to use some number to substitute for infinity.

inf:= 10: #Substitute for infinity.
Bcs:= theta(inf) = 0, D(theta)(0) = -k*(1-theta(0)):
ode:=
     diff(theta(eta),eta$2) +
     pr*Q*(1/alpha*(1-exp(-alpha*eta))*diff(theta(eta),eta) +
     Omega*theta(eta))+B*(-alpha*exp(-alpha*eta))^2 =
     0
;

Now you'll need two boundary conditions for the differential order of the ODE, plus one boundary condition for every unspecified constant: B, Q, Omega, k, pr, and alpha. So, you need eight boundary conditions.

Assuming that your matrix represents an "augmented" matrix, i.e., one where the rightmost column represents the right sides of the equations, then do

LinearAlgebra:-LinearSolve(R, free= t);

The t can be any unused variable that you want to use for the parameters. If you don't supply one, Maple will make one up.

You need to substitute a procedure (colloquially called a function) for x[1] rather than substituting an expression for x[1](whatever). Like this:

eval(junk, x[1]= ((a,b,c)-> R(b,c)*sin(omega*a + phi(b,c))));

Always use eval for these high-level calculus substitutions. Using subs isn't necessarily wrong, but eval is safer. The use of subs should be reserved for low-level purely syntactic substitutions.

I am highly skeptical of Kitonum's Answer. I don't doubt the computed critical value of D(theta)(0)---I got the same value. But I do doubt the solution curves. And I doubt the existence of a finite T such that theta(T) = Pi and D(theta)(T) = 0. My intuition is that diff(theta(t), t) (for t > 0) is a strictly positive, strictly decreasing function whose limit at infinity is 0.

Here are some computations to support these hypotheses:

restart:

I simplified your ODE. I hope that you'll agree that my form is equivalent.

ODE:= diff(theta(t),t$2) = sin(theta(t))*(cos(theta(t))-9.81);

diff(diff(theta(t), t), t) = sin(theta(t))*(cos(theta(t))-9.81)

The following dsolve solves the problem as a BVP. The value of the upper boundary T is up to you. As you make it higher, eventually either the BVP solver will give up, or you'll run out of patience waiting for it to finish.

T:= 9:

sol:= dsolve(
     {ODE, theta(0) = 2*Pi/3, theta(T) = Pi},
     numeric,
     abserr= 1e-9, maxmesh= 2^9
):

plots:-odeplot(
     sol, [[t,theta(t)-Pi], [t,diff(theta(t),t)]], 0..T,
     color= [red,blue], gridlines= false
);

Is theta'(T) = 0?

eval(diff(theta(t),t), sol(T));

HFloat(9.886613896682558e-13)

What is D(theta)(0)?

DTZ:= eval(diff(theta(t),t), sol(0));

HFloat(3.2496153618561823)

Now solve as an IVP using the high-accuracy taylorseries method.

sol2:= dsolve(
     {ODE, theta(0) = 2*Pi/3, D(theta)(0) = DTZ},
     numeric, method= taylorseries,
     abserr= 1e-9
):

My intuition is that the blue curve going negative simply represents the point where the Taylor series expansion is no longer accurate. If I ask for more accuracy, this point moves to the right.

plots:-odeplot(
     sol2, [[t,theta(t)-Pi], [t,diff(theta(t),t)]], 0..T,
     color= [red,blue], gridlines= false
);

 

 

 

Download Highly_unstable_ODE.mw

In general, there's no way to do it. But many builtin (or kernel) functions come from open-source libraries. For example, reading the code for isprime, one sees that it uses (for some inputs) gmp_isprime. That's part of the GMP open-source library. So, if you Google "gmp isprime", you'll quickly find that code in C.

How about saving everything to one .mla file?

save H||(1..n), "H.mla";

I just read your other Question, regarding Grid, and I now realize that you probably really need nine separate files. That can be done for the kth file by

save cat(H,k), sprintf("H%d.mla", k);

By the way, you may be confusing things by calling these .mla files. Use another extension. If you use .m, the files will be compressed into Maple's own text-based but not-human-readable format, which'll probably save an iota of time.

That needs to use nested seqs:

seq(seq([A[i],B[j]], i= 1..5), j= 1..4);

or

seq(seq([A[i],B[j]], j= 1..4), i= 1..5);

A better way, which doesn't require hardcoding the list sizes, is

seq(seq([i,j], j= B), i= A);

A method that works for an arbitrary list of lists L, each of arbitrary length, is

L:= [A,B]:
CP:= combinat:-cartprod(L):
'CP[nextvalue]()' $ mul(nops~(L));

 

You can't assign to a procedure's parameter, n in this case. You need to copy the parameter to a local variable, and then assign to that.  Like this:

cs:= proc(N::posint)::nonnegint;
description "Find the number of steps for a Collatz sequence to reach 1.";
local n::posint:= N, q::posint, r::identical(0,1), count::nonnegint;
     for count from 0 while n > 1 do
          q:= iquo(n, 2, 'r');
          n:= `if`(r=0, q, 3*n+1)
     end do;
     count
end proc:

Also, your count was off by one to my mind: cs(1) should be 0.

The properties (in the strict Maple sense of "properties") of an object can be inquired with the about command; however, this is only useful if you've used the assume command, which is the command which gives objects properties.

Somewhat similar to properties are attributes, which can be inquired with attributes. An object will only have attributes if you give it attributes with setattribute. You can give attributes whatever meaning you want.

To see all subexpressions of an expression e, use indets(e, anything). The anything can be replaced with any type expression (in the strict Maple sense of "type" (see ?type)). For example, all the integrals in e can be found with indets(e, specfunc({Int,int})). This makes indets one of the most useful and most commonly used commands.

To have e broken down into a logical (but not visual) tree which includes the smallest detail of information, use ToInert(e). This is often too much information to wade through by reading; however, this information can be manipulated programmatically. Some of this information may be arranged in a logical fashion that doesn't exactly correspond to the internal representation of the object, but does correspond to the logical representation of the object.

To have e broken down into a more-visual tree, but which uses some more-cryptic abbreviations, use dismantle(e). This information, as far as I know, corresponds precisely to the internal representation. This information is quite difficult to manipulate programmatically.

Ultimately, everything in a computer is stored as numbers. If you really want to go down to that level, use disassemble(addressof(e)). This will break it down one tree level only. Each given number can be broken down further with disassemble if it's an address. It's a good idea to follow the dismantle tree when you are doing this.

You should put your Maple version in the header of your Question. There's check-off boxes and a pull-down menu for it.

Your zero-one Matrix can be created like this:

n:= 4:
subsMatrix:= <seq(x[], x= Iterator:-CartesianProduct([0,1]$n))>;

Your list of substitution equations can be created directly from the iterator like this:

A:= [seq(a[i], i= 1..n)]:
subsEqs:= [seq(
     select(e-> rhs(e)=0, A=~ convert(x[], list)),
     x= Iterator:-CartesianProduct([0,1]$n)
)];

So note that there's no need to create the subsMatrix if its only purpose was to be an intermediary for the creation of the subsEqs.

Your error? It's probably that you need to say LinearAlgebra:-RowDimension rather than simply RowDimension and likewise for ColumnDimension.

The Iterator package comes with Maple 2016. If you have older Maple, let me know.

In bcs, you need to put parentheses around f1 in the first boundary condition and around f3 in the fourth boundary condition.

The command indets(sys2 union {bcs11} union bcs, name) can be used to find such problems. It this case, it told me that D was being considered a name (as well as an operator).

First 215 216 217 218 219 220 221 Last Page 217 of 395