Carl Love

Carl Love

26730 Reputation

25 Badges

11 years, 251 days
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

There is a way to solve this without resorting to commands outside of dsolve. The help page ?dsolve,numeric_bvp,advanced gives advice on how to apply dsolve's options to resolve some error messages given by `dsolve/numeric/bvp`. For messages related to Newton iteration convergence, it suggests the 'continuation' option, which is defined on help page ?dsolve,numeric,bvp (with two examples of its use given on the former help page).

To use the method, one introduces a parameter λ into the problem such that λ can vary continuously from 0 to 1, with λ=0 yielding a relatively easier BVP and λ=1 giving the desired BVP. The parameter can appear anywhere in the problem: in the ODEs, or the boundary conditions, or both. There is some subtlety required in deciding where and how to place λ. I tried this two ways with the problem at hand, and they both worked.

My first idea was to have λ=0 eliminate the convolution of f(x) and g(x), so I applied it to the g(x) term in the first ODE and to the term containing f(x) in the second ODE.

ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 ) (x) + lambda*g(x) = 0;
ode2:= (D@@2)(g)(x) + lambda*3*10*(f*D(g))(x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = 1, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color= black);

I recommend plotting the interval 5.7 .. 6 separately because there is an interesting fluctuation there which is not visible at the g-axis scale of the overall plot.

My second idea was to use λ=0 to "level" the boundary conditions so that everything starts and ends at 0 on the vertical axis. In this case, that simply means changing g(0)=1 to g(0)=λ.

ode1:= ( (D@@3)(f) + 3*f*(D@@2)(f) - 2*D(f)^2 + g ) (x) = 0;
ode2:= ( (D@@2)(g) + 3*10*f*D(g) ) (x) = 0;
bcs1:= D(f)(0) = 0, f(0) = 0, D(f)(6) = 0;
bcs2:= g(0) = lambda, g(6) = 0;
sys:= {bcs1, bcs2, ode1, ode2}:
dsn:= dsolve(sys, numeric, continuation= lambda);
plots:-odeplot(dsn, [x, g(x)], 0 .. 6, color = black);

Same plot, of course.

You are mixing new and old matrix commands. The new commands all begin with capital letters and tend to be spelled in full rather than abbreviated. Use Matrix and Vector instead of array. Use Adjoint instead of adjoint. Use Determinant instead of Det. See ?Matrix ?Vector ?Determinant etc.

I don't know why this unanswered question from seven years ago is appearing near the top of the "Active Conversations" list. Perhaps an answer was recently deleted and that counts as activity.

Anyway, I am guessing that you are trying to compute the exact determinant of a matrix over the integers. In that case you should use LinearAlgebra:-Determinant(A, method= integer). If your matrix is over some other ring, let me know.

After using the command kernelopts(opaquemodules= false), you can use the command showstat to list the code of any procedure in a module (or package), both exported and local procedures.

So, from your code, I'm guessing that for each integer i from 1 to 100, you want to find the smallest n such that tau(n) = i.  In other words, you want a partial inverse of the tau function. (Note to other readers: tau(n) is the number of positive integer divisors of n.) But your method won't work because for i = 97, n will have 29 decimal digits, and the whole age of the universe wouldn't be enough time to get to that number if you work by counting up from 1. I don't want to say what that 29-digit number is, because it would give too much of a hint.

Hint: You are correct that the prime factorization of n is the key to computing tau(n). Indeed, the primes themselves don't matter; it can be computed just from the exponents in the factorization. Take a look at showstat(numtheory:-tau). It's a very short procedure; only the last actual line of code matters here; the other lines are just there to quickly handle the trivial cases. Once you know what you want the exponents to be, just pick the prime bases that minimize the product. 

By adding sqrt(1+2*Pi), sqrt(1+4*Pi), sqrt(2*Pi), etc., to the bases used by identify (by using the extension option), I've identified several of the numerical solutions.  If cpu time becomes an issue, reduce the degree used by identify to 2 from its default 6. See ?identify

I'd like to take a look at this. What does it say when it crashes? What are the numbers for "Memory" and "Time" in the bottom edge, right side of your Maple window? What values do you get if you execute the commands kernelopts(cpulimit) and kernelopts(datalimit) in a fresh Maple session?

I took a look at your code, and I think that you're wasting cpu time by checking for small and medium factors with igcd, because isprime already does those checks. Take a look at showstat(isprime).

You set N:= (b-a)/h, with floating point a, b, h; so, N is floating point. But when you access N as the upper limit of a seq, it becomes integer. You should say N:= round((b-a)/h),

Just use eval(sol, _Z1= N) rather than eval(sol, _Z1~= N). A symbol with assumptions cannot be accessed by appending ~ to the name, even if you use name quotes (``); so eval(sol, `_Z1~`= N) will not work either. I don't think that the help is at all clear about this; it just shows an example of using indets to extract the name. 

You have made here a classic example of a programming mistake known as a "fencepost error." If a fence has 10 posts spaced one foot apart, how long is the fence? Nine feet, not 10. (Don't believe it? Then consider the case with just two posts.) Similarly a list with n elements has n-1 adjacent pairs of elements. The loop for i from 0 to c-1 do will execute c times.  If you change the lower bound from 0 to 1, then it executes c-1 times, which is the essence of Adri's brief answer.

If you continue with programming, memorize this bit of mental arithmetic, which you will find yourself using often in any system, not just Maple: The loop for i from a to b executes b-a+1 times, for integers a and b.

Could you state in more precise mathematical (set theoretic) language what you want?  Do you want one of these:

  1. all k in {$1..10} such that there exists a pair (a,b) in A^2 such that the relation holds,
  2. all pairs (a,b) in A^2 such that there exists k in {$1..10} such that the relation holds,
  3. all triples (a,b,k) in A^2 x {$1..10} such that the relation holds?

Also, as Doug pointed out, is the relation a^k = b or a^k = b^2? If the latter, then we can immediately eliminate all odd k.

I tried to download the files from your website, but I could not.  Maybe you could upload them directly. 

Does the fsolve work at, say, Digits = 15?  Try using the fulldigits option to fsolve at your higher value of Digits, 90.


This behavior is not specific to matrices. The object created by making assignments to a subscripted (a.k.a. indexed) name is always a "table". If a list (or Vector, Array, etc.) already exists, then the entries can be accessed by subscripting.  The most common way to create a list is to create a sequence and enclose it in square brackets ("[ ]"). For your example:

M:= [seq(Matrix(2,2, [[1, i], [i, i^2]]), i= 1..4) ]

Function application (f(...)) automatically (i.e. inherently to the syntax) distributes over list (or set ({})) creation, so the above could be

M:= [seq] (Matrix(2,2, [[1, i], [i, i^2]]), i= 1..4)

which I much prefer the look of, especially when the seq command spans several lines of code.


One may be tempted to create a list by repetively appending to an existing list, as in

M:= [];

for i to 4 do

    M:= [M[], Matrix(...)]


This does create M as a list, but the problem is that it creates and discards 4 temporary copies of M.  The time complexity for this will be O(n2) where n is the list length. This will be a problem if n is large.  If one must use a for loop to create the list's entries, it is best to store them in a table within the loop, then convert the table to list outside the loop:

for i to 4 do  M[i]:= Matrix(...)  od;

M:= [seq](M[i], i= 1..4);





[1,2] + [3,4] *~ 3 *~ (x-1)^2;

Enter your integral in unevaluated form, i.e. as Int with a capital "I", and with symbolic e_0 as the upper limit:

t:= Int((1+121/304*e^2)^(1181/2299)*(1-e^2)^(-3/2)*e^(29/19), e= 0 .. e_0);

Then, simply,

plot(t, e_0= 0 .. 0.9);


Note to other readers: In this case, I needed to make sure that Maple did not try to evaluate the integral symbolically.  There is a very serious bug such that evaluating the integral (i.e. value(t)) causes the Maple kernel to utterly and completely crash after using about 11 seconds of processor time.

First 382 383 384 385 Page 384 of 385