Carl Love

Carl Love

28085 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

The following simple code is probably adequate for your needs. It could be made slightly more efficient by pre-sieving for the primes.

MyProperty:= proc(q::prime, P::name) 
local p:= ithprime(q - numtheory:-pi(q));
    if p >= q then P:= p; true else false fi
end proc
:    
MySeq:= proc(n::posint)
local k:= 0, R:= Vector(n), p, q:= 1;
    while k < n do 
        q:= nextprime(q);
        if MyProperty(q, 'p') then
            k:= k+1;
            R[k]:= [p,q];
        fi
    od;
    [seq](R)
end proc
:
S:= MySeq(9);
S := [[2, 2], [13, 11], [17, 13], [29, 17], [31, 19], [43, 23], 
  [67, 29], [71, 31], [97, 37], [107, 41]]

(p_seq, q_seq):= map(op~, [1,2], S)[];
   p_seq, q_seq := [2, 13, 17, 29, 31, 43, 67, 71, 97, 107], 
     [2, 11, 13, 17, 19, 23, 29, 31, 37, 41]


 

In Maple, type checking can only be done at run time. This is mostly because types can be arbitrarily complex, and checking them may involve execution of an arbitrary amount of Maple code. Of course, for the example that you show, it'd be theoretically possible before runtime; but I'm sure that you realize that your example is extremely superficial. 

The code below allows you to set the values of the parameters with sliders and have the plot update automatically. For some values of the parameters, there will be warning messages such as "cannot evaluate the solution further (left/right) of ... probably a singularity." These warnings do not prevent you from proceding with the Explore. If you want, the field arrows could be added.

restart:.
A:= a_1*c^3 - b_1:  B:= 3*a_1*c^2 - 3*b_1/c:  C:= d*c - e:
sys:= {diff(x(t), t$2) = (B*x(t)^2+C*x(t)-e)/A, x(0)=f, D(x)(0)=g}:
P:= [a_1, c, b_1, d, e, f, g]:
sol:= dsolve(sys, numeric, parameters= P):
MyPlot:= proc(P)
    sol(parameters= P);
    plots:-odeplot(
        sol,
        [x(t),diff(x(t),t)],
        t= -3..3,
        view= [-6..6, -6..6],
        thickness= 2, axes= boxed, color= red, scaling= constrained
    )
end proc
:
Explore(MyPlot(P), parameters= (P=~ -5..5));

 

subsindets(B, string, parse)

As Tom admits, explicitly entering the expression to be integrated is "ugly", and I find Axel's modification to be only a slight improvement. And that ugliness is made worse by the fact that the OP actually has a large number of such sets of integrals to do that I removed from the posted worksheet because the syntax issue is identical even though the integrands are not.

I like to build the expression to be integrated directly into the ODE system, like this:

restart;
eq1:=diff(f(y), y$4)+Uhs*diff(E(y),y$3)-(diff(f(y), y$2))+(diff(theta(y), y$1))= 0:
eq2:=diff(theta(y), y$2)+(diff(f(y), y$2)+1)^2+1+diff(theta(y),y$2) = 0:
E:=y->zeta*(cosh(k/2*(h1+h2-2*y)))/(cosh(k/2*(h1-h2))):
bcs:=f(h1) = -(1/2)*(Q-1-d), f(h2) = (1/2)*(Q-1-d), (D(f))(h1) = -1, (D(f))(h2) = -1,theta(h1) = 0, theta(h2) = 1:
epsilon1:=0.1:d:=1:omega:=Pi/6:h1:=-(1+epsilon1*sin(2*Pi*x)):h2:=d+epsilon2*sin(2*Pi*x+omega):F:= Q-1-d:epsilon2:=0.5:x:=1:
alpha:=Pi/6:
de:=eq1,eq2,bcs:
d1 := subs(Uhs =-2, zeta=3,k=1,{de}):
param:= {Uhs =-2, zeta=3,k=1}:
P1:= eval(diff(f(y), y$3)+Uhs*diff(E(y),y$2)-(diff(f(y), y$1)+1)+(theta(y))+sin(alpha),param):
ec:=0.5:
for Q from -3 to 3 by ec do 
    F2[Q]:= dsolve(
        d1 union {diff(J(y),y)=P1, J(h1)=0}, 
        numeric, maxmesh= 25500, continuation= lambda1
    ) 
end do:
for Q from -3 to 3 by ec do 
    P3[Q]:= eval(J(y), F2[Q](1)) - eval(J(y), F2[Q](0))  
end do;
                   P3[-3] := 2.30876985489011
                  P3[-2.5] := 1.36946865887356
                 P3[-2.0] := 0.460748411800256
                 P3[-1.5] := -0.417186862674401
                 P3[-1.0] := -1.26412891537822
                 P3[-0.5] := -2.07986512872413
                  P3[0.] := -2.86417840284371
                  P3[0.5] := -3.61684700714295
                  P3[1.0] := -4.33764445144679
                  P3[1.5] := -5.02633933429334
                  P3[2.0] := -5.68269517762163
                  P3[2.5] := -6.30647027769476
                  P3[3.0] := -6.89741755441151

 

You can place the cursor wherever you want (such as after the x in dx), press backspace, and change it to y. Likewise, place the cursor before the "d" to enter/edit the integrand, then start a new definite integral there. There's no need to use parentheses, but I suppose that you could if you really wanted to. In the hyperlink posted below, I have entered the expression which in standard Maple Notation would be

Int(Int(f(x,y), x= c(y)..d(y)), y= a..b)

https://learn.maplesoft.com/index.html#/?d=DQJKHPFIKPOTBIGPKTGHCFFPPTKSDPBSDGLKFNFUFIGUKNBQAFDLJOOKNHDUERMQKKCSJKPGASGFEIPQLNCUIMMRHKBFMHHGANER

I think that the problem is simply that you've made a mental arithmetic error in determing your expected result, because x^(2*k)/x^k is x^k, not x^2.

An unrelated tip: You shouldn't assume that the bound variable of the Sum is literally (even if that's always true). Instead, extract it via op([2,1], s) or, equivalently, lhs(op(2,s)).

The code that acer wrote for you that uses operators such as #mi(...) and #mo(...was intended to display output via Maple's GUI, not via Latex. I vaguely recall him explicitly saying that, but upon rereading, I can't find it, so maybe he didn't say it. Anyway, I'm certain that that was his intention. Those operators are part of the language HTML (I believe), and not the language Latex.

Now I wouldn't be at all surprised if Latex had some facility for interpretting HTML. If it does, it's likely possible to write a brief Maple procedure to edit Maple's Latex output in such a way that that facility was turned on.

From my personal observation of the expansion of text-formatting languages, there was a rapid divergence about 25 years ago between those primarily intended for paper documents (such as Latex and PDF) and those primarily intended for on-screen documents (such as HTML and XML). Note that I said "primarily": Almost all of them can be used on either medium, but for almost all of them there's one medium that's distinctly much easier. I'd be quite interested in reading here anyone else's opinion on this divergence.

Your main error is attempting to define a function (or procedure) using this syntax:

f(...):= ...;  #bad way

Unfortunately, the above is not an outright error, and Maple allows you to "get away with" doing this to some extent, and then the new user becomes habituated to it. But it doesn't always work, and, to my mind, should never be used to define a procedure (the above syntax has other valid and perfectly acceptable uses that aren't further discussed in this article).

There are several correct ways to define a procedure:

  1. f:= x-> e #univariate case
    where e is any explicitly specified expression containing x.
  2. f:= (x, y, ...)-> e 
    allows any number (including 0) of parameters (aka variables).
    must be still be an explicitly specified expression of the parameters.
  3. f:= proc(x, y, ...) s1; s2; ...; e end proc
    This allows for arbitrary length and number of instructions. e doesn't quite need to be explicit, but expressions created outside the procedure using the procedure's parameters are not allowed.
  4. f:= unapply(e, [x, y, ...])
    Now e doesn't need to be explicit, and it can be created in a separate place.
  5. f:= subs([x= e1, y= e2, ...], eval(P))
    where is an existing procedure containing undeclared global variables x, y, ..., and e1, e2, ... are any expressions with almost no restrictions (they can even use P's parameters).

And there are even further ways that I'd classify as "meta-programming".

So, here are 4 ways to perform your task, in increasing order of sophistication. Actually, the individual components of each of these are fairly independent, so you can mix-and-match.

Download func_norms.mw

MaplePrimes won't expand my worksheet inline, but you should download it from the link immediately above. Nonetheless, here it is with some of its output reformatted, and the plots removed (simply for brevity; anyway, they're not interesting plots; the interest lies in the code that creates them).

#Method 1: Doesn't create procedures at all.
restart: 
f:= x^2:
n:= int(f^2, x= 0..1):
f:= f/n;
                                          2
                                  f := 5 x 

plot(f, x= 0..1);

#Method 2: Creates procedures with -> and unapply. 
restart:
f:= x-> x^2:
n:= int(f^2, 0..1):
f:= unapply(f(x)/n, x);
            f := x -> 5*x^2

plot(f, 0..1);

#Method 3: Uses a proc procedure to abstract (slightly) the concept of 
#"functional normalization". This procedure optionally displays the plot
#as a side effect with the "normalized" function as the return value.
#
#Attn nitpickers: Yes, I know that this isn't a true norm. I'm just
#working with the OP's formulation. A true norm is covered in Method 4.
restart:
normalize1:= proc(f, ab::range(algebraic), N)
local x, r:= unapply(f(x)/int(N@f, ab), x);
    userinfo(1, procname, NoName, print(plot(r, ab)));
    eval(r)
end proc
:
#Set next line to 1 if you want to see the plot, 0 if you don't:
infolevel[normalize1]:= 1:

f:= x-> x^2:
f:= normalize1(f, 0..1, f-> f^2);
            f := x -> 5*x^2 

#Method 4: Further abstraction of the normalization concept so that the
#whole norm can be specified as a procedure.
restart:
normalize2:= proc(f, N, {Popts::list:= []})
local x, r:= unapply(f(x)/N(f, _rest), x);
    userinfo(1, procname, NoName, print(plot(r, Popts[])));
    eval(r)
end proc
:
infolevel[normalize2]:= 1:
f:= x-> x^2:

#Create a standard p-norm (only standard if p>=1):
N:= (f, ab::range(algebraic), p::algebraic)-> 
    int((f-> f^p)@abs@f, ab, 'continuous')^(1/p)
:
f:= normalize2(f, N, 0..1, 2, Popts= [0..1]);
         f := x -> x^2*5^(1/2) 

infolevel[normalize2]:= 0: #Turn off plotting.
f:= x-> x^2:
f:= normalize2(f, N, 0..1, p); #p is now symbolic.
    f := x-> x^2/(1/(2*p+1))^(1/p)

#Instantiate p (i.e., give it a specific value):
eval(f(x), p= 2);
                                   2  (1/2)
                                  x  5     

#Another way to instantiate: subs into an existing procedure:
#This is not that useful in this particular example. I just wanted to
#briefly illustrate the 5th method of procedure creation that I itemized at the 
#beginning of this article.
f:= subs(p= 2, eval(f));
         f := x -> x^2*5^(1/2)

#But be aware that the "subs" method is quite tricky to master, yet 
#extremely powerful once you do so.


 

See ?readlib. I think that this command has been obsolete for well over 20 years. But older usages, such as some appearing in dsolve, are still in some library code. Calls to readlib do nothing, and they're not needed anyway, so they're harmless.

This just requires a simple wrapper that resets assertlevel twice with the call to FetchAll in between:

MyFetchAll:= proc()
local 
    al:= kernelopts(':-assertlevel'= 1), 
    r:= Database:-SQLite:-FetchAll(args)
;
    kernelopts(':-assertlevel'= al);
    r
end proc
:    

Boundary-layer ODE BVPs like this related to nano-fluid dynamics are one of the most common sources of Questions here on MaplePrimes (perhaps even the single largest source), and I've worked through several such papers here. In particular, see my Post "Numerically solving BVPs that have many parameters" . Here is how to read this paper with an eye towards computing the solutions in Maple. The parts of the paper that you actually need to use are extremely short.

1. Everything in the first two sections (10 pages) of the paper can be ignored, except to the extent that after the solution procedure is coded, we'll need to look in those sections for the ranges of numeric values of the parameters to use.

2. The complete BVP such as it's needed for computation is given in Eqs. 21 - 23 on page 11. That's all there is! Here, the authors have renamed the dependent variables of the system as y(1), ..., y(6)yy__1yy__2. That's a complete load of crap! Just stick with the original names f(eta), ..., diff(f(eta), eta$4)theta(eta), ..., diff(theta(eta), eta$2), which are the names commonly used in nano-fluid dynamics (as I'm sure you're aware, this being not the first time that you've posted to MaplePrimes about such a BVP). If the authors have done this renaming because Matlab syntax requires it, then shame, shame on Matlab!, because this notation is hideous. Maple handles this conversion to a first-order system in the background.

All parameters that appear symbolically in Eqs. 21 - 23 must remain symbolic in your coding of the BVP at this point. That is, we won't assign them numeric values until dsolve is called.

3. One of the boundary conditions is stated "as eta -> infinity". Maple's BVP solver doesn't allow for that. We simply pick a value (such as 1) to use instead of infinity. But don't directly use 1! That would be a very bad programming decision, likely to lead to confusion, error, and/or inflexibility. Instead, give this "fake infinity" a name (such as inf), and use that exclusively for the rest of the code. Thus, that 6th boundary condition should be theta(inf) = 0. It'll be worthwhile to look through the paper more thoroughly to see how the authors handled this. 

4. Everything up until this point is extremely straightforward, and the details have been shown hundreds of times (at least) on MaplePrimes.

5. The major hurdle with these boundary-layer BVPs (and the cause of a great number of those MaplePrimes Questions) is a failure to converge for some values of the parameters due to non-uniqueness of the solution. (In my limited experience, this is more likely to happen when a "fake infinity" is used.) A method to handle this is discussed in the first paragraph of section 4. The syntax for doing this in Maple uses the approxsoln option to dsolve.

It will be interesting to see if Maple can perform as well as, or better than, Matlab for this. From a computational viewpoint, this is the only interesting part of the problem, because, as I said, everything else is completely straightforward. But this final part is quite interesting.

The entirety of this process can be done without the slightest bit of knowledge of fluid dynamics; indeed, I have no such knowledge myself. It's a purely a Numerical Analysis problem.

@Wilks Yes, what you just said about restart is correct, and thus you can't directly use Kitonum's Answer above. I think that he was assuming that Y1 and its derivative were undefined functions.

Assuming that your `Y1"`(x) is a correctly defined function (or algebraic expression) for which Maple can find a symbolic antiderivative, your C1 can be expressed simply as 

-intat(`Y1"`(x), x= 0)

If you prefer, you could also write explicitly C1 = -intat(`Y1"`(x), x= 0). Regardless of the complexity, there'd never be a need to use solve for this; however, if for some reason (pedagogy, perhaps) you'd prefer to use solve, you could do

solve({intat(`Y1"`(x), x= 0) + C1}, C1)

If you want a procedure to do this for arbitrary initial conditions, you could use this:

IC:= proc(`y'`::appliable, ic::function(algebraic)=algebraic, C::name)
local x;
    `tools/genglobal`(`if`(nargs=3, C, ':-_C')) = 
        rhs(ic) - intat(`y'`(x), x= op(lhs(ic)))
end proc:

And you'd use it like this (assuming that `Y1"` has been defined, Y1 has not been defined,  x0 and y0 are any appropriate algebraic expressions (including simply numbers), and C1 is an unassigned name):

IC(`Y1"`, Y1(x0)=y0, C1);

where the 3rd argument, C1, is optional, and if not supplied the procedure will choose an unused name beginning with _C.

The command forget can be used to dump the cache tables and remember tables. See its help page ?forget. However, in my opinion (based on many years experience intensively time testing Maple code), more meaningful results can be obtained by calling the procedure multiple times with random input. Generate a random list of sets of polynomials (of the type that you're interested in), and map the procedure over that list.

Use CodeTools:-Usage, and note both the real time and cputime. 

Note that your desired unexpanded form f^6/6 has a constant term of 2048/3. This constant term doesn't appear in Maple's int result, so that result can't be factored to f^6/6.

First 80 81 82 83 84 85 86 Last Page 82 of 395