Carl Love

Carl Love

21198 Reputation

24 Badges

8 years, 246 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

Yes, you are correct that your piecewise expression can be converted to an equivalent form that doesn't use piecewise and this makes some improvement to the time of the integration. After defining the piecewise via f:= piecewise(...), do

f:= convert(f, Heaviside):

Other recommendations:

Change Digits:= 30 to Digits:= 15, which is the largest value for which hardware-float computations will be used, which are much faster. As far as I can see, this doesn't change the final result, but my observations of this aspect have been superficial thus far.

Remove Digits:= 4. In my opinion, that's too low to ever be used as the global value of Digits. There are other ways to control the precision of a numeric integration, which are discussed in the next paragraph.

Inside the int command, after keyword numeric, add the options epsilon= 0.5e-4, digits= 7. The returned result will be guaranteed to 4 significant digits. There is some chance that the integral will be returned unevaluated, which means that it wasn't able to control the accuracy sufficiently. If this happens, let me know.

If you make these changes, then I think that the majority of the computation time is consumed by solve, not by int. If this amount of time is still too much for you, let me know.

The sequence of iterates in Newton's method can be extremely sensitive to the number of digits of precision at which the calculations are performed. This is totally normal.

You may already know this, but that is not conveyed in your Question or worksheet. Since the second derivative of exp(x) never changes sign (in other words, there are no inflection points), there can be at most two intersections between the curve and any straight line, and if an intersection is also a point of tangency, then it's necessarily the only intersection. The cases of 0 intersections and 1 intersection are relatively easy to classify, as has already been done. Thus, all remaining cases are 2 intersections.

Do you mean this?:

data:= convert(
    Statistics:-Sample(Uniform(-2, 2), [100, 2]), 
    list, nested
    color= COLOR~(
        .5 -~ argument~(-(Complex@op)~(data)^~2)/~(2*Pi)
    style= point, symbolsize= 20, symbol= solidbox

Maple has two packages for parallelization: Threads and Grid:

  • With Threads, the procedures running on different cores share all memory except for variables that you explicitly tell them not to share. 
  • With Grid, the procedures in different cores run as distinct "processes", each with its own memory, and they can only communicate with each other (if any communication is needed at all) via variables or messages that you explicitly send between them.

The vast majority of Maple's high-level symbolic library commands make extensive use of global variables, which makes them unsuitable for the shared-memory environment of Threads. I suspect that Query is in this category.

Try using the command Grid:-Map for this.

I'm guessing that you might want some tutorial help with doing this problem "by hand".

The first step is to "complete the square" in the denominator, thus obtaining 1/((s+4)^2 + 4). We recognize this as the composition f@g of the two functions g:= s-> s+4 and f:= s-> 1/(s^2+2^2). The inner function g is a "shift": just adding a constant a to s. That means the inverse transform will contain exp(-a*t) as a factor (not as a function composition). The remaining part is of the form 1/(s^2+b^2). From a table of transforms (even the most-basic tables will likely have this one), the inverse transform of that is sin(b*t)/b. This is the other factor, making the answer exp(-4*t)*sin(2*t)/2.

A second approach to the problem is much more straightforward---it doesn't require much recognizing of things as being of certain special forms---but it does delve into complex numbers. We factor the denominator over the complex numbers to obtain (s+(4+2*I))*(s+(4-2*I)). Compute the partial fraction decomposition of the rational function from that factorization to obtain 

I/4*(1/(s + (4+2*I)) - 1/(s + (4-2*I))

Each fraction is a "shift", so the inverse transform is

I/4*(exp(-(4+2*I)*t) - exp(-(4-2*I)*t))

If you convert that to real form (which can be done with Maple command evalc) you get the same thing as by the previous method.

Your transcription of your code is so chock full of typos that the only reason that I can understand your Question is because I worked on an earlier Question of yours to which this one is closely related. I suspect that very few, perhaps 0, other readers can understand it at all, which may be why it's taken a long time for you to get an Answer.

In addition to the typos, the expressions -Pi/2 + 10 and Pi/2 - 10 are surely wrong. By 10, did you perhaps mean 10 degrees? In that case, you need to change 10 to Pi/18.

Now let's address your undefined issue: As is well known and explicitly taught---even in first-semester calculus---the derivatives of piecewise functions are often undefined at points where one piece changes to the next. So the fact that you're getting undefined in your piecewises is not an error or a bug (and I'm not saying that you were claiming that it was an error or a bug). But I understand that some further symbolic processing that you're doing with these expressions is having problems because of the undefineds. That's totally believable. The following example serves to illustrate where these undefineds come from, and also shows a method to avoid them. Please execute the code:

f1:= abs(x):
f2:= convert(f1, piecewise, x);
diff(f2, x);
#Note the correct presence of undefined in previous result.

f3:= convert(f2, Heaviside);
convert(diff(f3, x), piecewise, x);
#Same result as before

Heaviside(0):= 0:
convert(diff(f3, x), piecewise, x);
#Now there's no undefined.

There are many subtleties involved in working symbolically with piecewise expressions, especially when there is more than one independent variable (it looks like you have both x and y), so I'm not claiming that the above will be a complete solution to the undefined issue. See if you can make *some* progress with it, and please report back.

I agree that the result 1 doesn't make sense and should be corrected. But compare your result with this:

convert(abs(1,x), Heaviside) assuming x::real;
                      -1 + 2 Heaviside(x)

This small example, as well as Tom's help page excerpt, should answer your Question about $:

f1:= proc(x, y) x^2+y^2 end proc:
f1(a, b, c);
                             2    2
                            a  + b 

f2:= proc(x, y, $) x^2+y^2 end proc:
f2(a, b, c);
Error, invalid input: too many and/or wrong type of
arguments passed to f2; first unused argument is c

The x and y in those examples are called parameters, not arguments. The ab, and c are the arguments.

Do not think in terms of grids or cells! Coloring this plot is a purely algebraic problem, and I can completely describe it at the following high level of abstraction so that it'll work for any two (suitably smooth and close-enough-to-homeomorphic) parameterized surfaces, in any two coordinate systems, on any two distinct grids (perhaps even of different size). The beauty of plot3d's syntax (especially the part related to color functions) is that it handles those grid-and-cell-level details for you in the background so that you can focus on the algebra.

Only the algebraic details matter here, so I'm no longer going to mention any topological (e.g., "homeomorphic") or analytic ("smooth") details. The algebraic details that matter most are the number of input and output dimensions of all of the following functions and the order in which those dimensions are given: the coordinate transformations, the parameterized surfaces, the projections, and the coloring functions. 

Define a 3D coordinate system as a surjective function C: R-> R3 such that C(a, b, c) =  <x(a,b,c), y(a,b,c), z(a,b,c)>, where <x, y, z> are the usual cartesian coordinates and <a, b, c> are the new coordinates. (In practice, most useful coordinate systems are named, predefined, and known to Maple, so we don't need to actually supply the function C.)  Define a parameterized surface S in coordinates C as a function S: A -> R3 such that S(u,t) = <a(u,t), b(u,t), c(u,t)>, where A is some 2D subset of R2, <a, b, c> are the coordinates of C, and <u, t> are the parameters of S.

For the time being at least, I'll only use 1 color system; let's say HSV, because that's what I used for the final example. Like most color systems, it has 3 coordinates, which has nothing to do with geometry; rather, it's based on the biology of (the most-common version of) human retinae; other animals often have a different number of dimensions of color vision. Naturally, I'll call the coordinates of HSV <h, s, v> (hue, saturation, value). For the purpose of use with plot3d, all 3 are bounded in the interval 0..1[*2]. Define a coloring function F for surface S as F: A -> (0..1)3 such that F(u,t)= <h(u,t), s(u,t), v(u,t)>, where A, u, t are as they were defined for S.

Given two surfaces S1 defined on A1 and S2 defined on A2, define a projection P from S2 to S1 as a (sufficiently bijective) function from A2 to A1[*3]. Note very carefully that the projections that I'm using here have 2 (not 3) input and output dimensions because surfaces are fundamentally 2D objects embedded in 3D space. Now let's say that we have S1, S2, and P (from S2 to S1), and we've decided on a coloring function F for S1. The algebraic trick that this entire article has been leading up to is that The coloring function to use for S2 is the functional composition F@P. It's as simple as that!

Now, onto your sphere-on-the-plane example. I chose to represent the sphere in spherical coordinates and the plane in cylindrical coordinates because that makes the algebra for all the functions easy enough to do in my head. For the sake of having 6 distinct[*1] coordinate variables, I represent spherical coordinates as <rho, psi, phi> where rho is the radius, psi is the longitude (in range -Pi..Pi), and phi is the angle between the radial vector and the positive z-axis (in range 0..Pi); and cylindrical coordinates as <r, theta, z> (which I probably don't need to define for you). The surface functions and the projection from the plane to the sphere are

Sph:= (psi, phi)-> <1, psi, phi>:  
Pl:= (r, theta)-> <r, theta, -1>:
P__pl_to_sph:= (r, theta)-> (theta, Pi-2*arctan(r)):

In other words, the projection is simply psi = theta, phi = Pi - 2*arctan(r). [Edit: I now know that this is not the projection that you intended; however it still works as an example here, and my algebraic/functional formulation allows this one-line detail to be changed without changing anything else.] Note that it's essential to be consistent about the orders of the 2 parameters (..., ...) on the left sides of -> and the orders of the 3 coordinates <..., ..., ...> on the right sides of ->, and this is usually (for me at least) the most difficult part of working with alternate coordinate systems. (However, the vectors <...> could be changed to lists [...]; that makes no difference to (modern versions of) plot3d. Also, Maple's main plotting commands don't care what variable names you use for the coordinates in any coordinate system.) Furthermore, the order of the coordinates for any of Maple's predefined coordinate systems must be obeyed, although it's poorly documented. 

For the coloring function of the sphere, I'm using

Color__sph:= (psi, phi)->
    COLOR(HSV, (psi/Pi+1)/2., 1.-phi/Pi, 1.-phi/Pi)  #note the decimal points

I expect that you'll make that more elaborate, but it's the only one of these functions that you should change. [Edit: See the Edit comment in the previous paragraph.] Finally, the plotting command is

        [Sph(psi, phi), Pl(theta, r)], 
        [psi, theta]=~ -Pi..Pi, [phi= 0..Pi, r= 0..4],
        coords=~ [spherical, cylindrical],
        color=~ [
            Color__sph(psi,  phi),
            (Color__sph @ P__pl_to_sph)(r, theta)
    lightmodel= none, scaling= constrained, axes= none,
    viewpoint= circleleft

I've tried several times to upload the resulting animation as a GIF file, but MaplePrimes is having a problem with that.

[*1] It's not actually necessary that the coordinate names between the coordinate systems be distinct. Making them distinct simply helps with my mental organization of the problem.

[*2] Bounding the color coordinates in the interval 0..1 is not strictly required. Maple will usually normalize the values for you. However, I don't like to rely on that; sometimes I've had problems if I don't normalize them myself.

[*3] This is perhaps not a standard definition of projection, and it's certainly not how that word is used in linear algebra, but it does conform both to the way that you've already used the term (as in "stereographic projection") and to how it's used in relation to geographic maps. 

A good solution is to create a Maple table-based index mapping the first two columns (which'll henceforth be called the key fields) to the row numbers (rows are henceforth called records). This solution is general, efficient, and remarkably simple---the comments for the code below are far longer than the code itself. In using this method, you'll never need to think about, see, or specify any record numbers; that's handled completely by the code. I know that you literally said that you wanted "to print the row numbers"; however, after you see this method you may realize that that's a totally unnecessary intermediate step rather than your final goal. Anyway, the code below will still easily let you see the row numbers if you really want to. Just do 


where DB is the return value of the procedure Index, below, and comb is exactly like you had in your Question.

Disambiguation of index and its inflected forms: In database theory, the usage of index and its inflected forms is more akin to normal English (as in, "the index of a book") than it is in general computer science and mathematics ("an index of a matrix or array"). Throughout this article, the noun index refers to the entire table-stored mapping between keys and record numbers; whereas in Maple (and other computer languages) it usually refers to a specific key value, as in the Maple command indices. The verbs index and reindex refers to the creation of an index rather than storing or looking up a specific key value (as in the Maple command index), for which the verbs store or look up will be used. The past participle indexed can be applied either to the entirety or to a specific case. The common-English verb key, its past participle keyed, and the adjective key are, to my knowledge, not used in any ways specific to databases, mathematics, or computer science; only the noun key is used.

Importing your data: I'll assume that you can import your data file into Maple as a Matrix. If you need help with that aspect, then ask. If the file is too big to fit into main memory as a Matrix, let me know. As long as just the index can be somehow squeezed into memory, a modification of this code can still work.

A matrix-based database implementation
We're given a Maple Matrix M each row of which represents a record of data. The columns are the fields. Some fields are key fields---those by which the records can be specified, found, or indexed. The key fields have a specified order, which needn't be their order as columns. So, if the key fields are represented by a list K of column numbers, then the key for the record at row r is M[r, K].  (The italicized words in this and the preceding paragraphs are usages specific to database theory.)

This implementation allows for a key to specify multiple records. The code below provides the option (multikey) of creating the index such that it'll retrieve all records that have a given key. Using this option incurs a substantial performance penalty when the index is created, although I think it's still fast. There is no penalty for the retrievals.

This implementation is designed to be efficient---regardless of the number of records---as long as adding new records or changing keys are relatively rare compared to retrieving records. There's no performance penalty for changing the data in non-key fields. Record deletion can be done efficiently by marking certain non-key fields as invalid data. Adding records or changing keys requires reindexing with Index.

Create an index for a database represented as a matrix M with each 
row representing a record.  K is a list of the key columns. 

The return value is a Maple "Record"---essentially an ordered pair 
[Data, Index] where Data is the matrix M and Index is the index
created by this procedure.  Such an ordered pair will be
called a "Database".  Note that "Record" is a stock Maple data
structure, and in this case it doesn't mean "record" in the sense
of database theory as used in the introduction.

Keyword option 'multikey':
Keys are allowed to point to multiple records.  If 'multikey' is
specified, then they will be indexed accordingly.  Otherwise, only
the last record for each key is indexed.  This option incurs a 
substantial performance penalty during the indexing, but no penalty
for the retrievals.

(*  This procedure creates a so-called 'sparse' table as the primary
data structure for the index.  A 'sparse' table simply means one 
that returns 0 instead of an unassigned reference when it's passed
a key that hasn't been stored.  You'll see two references to these
0s in the following procedure, 'Get'.  This works as an "invalid 
key" flag because 0 is never a valid row number.  *)

Index:= proc(
    M::Matrix, K::list(posint):= [1], {multikey::truefalse:= false}
local i, J, n:= upperbound(M)[1]; #n is # of rows
    if multikey then
        J:= table(':-sparse');
        for i to n do J[seq(M[i,K])][i]:= () od;
        J:= indices~(J, ':-nolist')       
        table(':-sparse', [seq](seq(M[i,K])= i, i= 1..n))
    Record("data"= M, "index"= (()-> J[args]))
end proc
Get the records (rows) in database D corresponding to the keys in K.
Each entry in K must be a "tuple" (i.e., a list) [key1, key2, ...].

The return value is the submatrix of the data corresponding to the keys.
If K is a list, then the rows are ordered accordingly.  If K is a set,
then natural order is used.

Keyword option 'verify':
If the optional keyword 'verify' is given, then if any of the keys 
are not found then those that are not found are returned.  Otherwise,
the return is the same as usual.  If there's likely to be unfound
keys, and 'verify' is specified, then it's slightly more efficient to 
make K a set.  
Get:= proc(
    D::record(data::Matrix, index::procedure),
    {verify::truefalse:= false}
local R:= (D:-index@op)~(K); #rows corresponding to the keys K
    if verify then
        `if`(0 in R, select(k-> D:-index(k[])=0, {K[]}), D:-data[[R[]]])
        D:-data[subs(0= (), [R[]])] #The subs removes 0s, if any
end proc

r:= 2^16: #number of rows
c:= 5: #number of non-key data columns
(k1,k2):= (1..2^9, 1..2^9): #ranges for keys 1 & 2

K:= LinearAlgebra:-RandomMatrix((r,2), generator= rand(k1)):
NK:= LinearAlgebra:-RandomMatrix((r,c), generator= rand(-2. .. 2.)):
M:= <K | NK>:

DB1:= CodeTools:-Usage(Index(M, [1,2], multikey)):
memory used=172.38MiB, alloc change=238.00MiB, 
cpu time=2.19s, real time=1.47s, gc time=1.86s

DB2:= CodeTools:-Usage(Index(M, [1,2])):
memory used=28.84MiB, alloc change=-8.00MiB, 
cpu time=468.00ms, real time=393.00ms, gc time=312.50ms

(* Creating multiple indices for the data does not duplicate the data. All
data remains stored in, and only in, its original matrix M, and it can be
edited in M. *)

#Create some keys to lookup:
combos:= convert(
    LinearAlgebra:-RandomMatrix((999,2), generator= rand(k1)), 
    list, nested
CodeTools:-Usage(Get(combos, DB1, verify)):
nops(%); #count invalid keys.
memory used=4.74MiB, alloc change=0 bytes, 
cpu time=15.00ms, real time=15.00ms, gc time=0ns


CodeTools:-Usage(Get(combos, DB1));
memory used=4.60MiB, alloc change=0 bytes, 
cpu time=0ns, real time=14.00ms, gc time=0ns
[240-row matrix output omitted]

CodeTools:-Usage(Get(combos, DB2));
memory used=4.59MiB, alloc change=0 bytes, 
cpu time=16.00ms, real time=13.00ms, gc time=0ns
[213-row matrix output omitted]




If you read through the code showstat(ifactor), you'll see that the specification of a method in the second argument is open ended. In other words, any name meth can be given, any then ifactor will pass your number on to `ifactor/meth`. So, you're not restricted to the methods mentioned at ?ifactor. As I said earlier, the Maple 2021 library has 30 procedures so named, and the tags meth can be extracted like this:

Meths:= cat~(``, index~(LibraryTools:-PrefixMatch(`ifactor/`), 9..-3));
[PollardP1, PollardP1/ext, PollardP1/tree, PollardRhoBrent, 
  PollardRhoBrent/ext, QuadraticSieve, SmallFactors, 
  SmallFactors/ext, SmallFactors/lib, anm1, easy, easyfunc, 
  from_signature, ifact0th, ifact1st, ifact235, lenstra, mixed, 
  morrbril, mpqs, mpqsmixed, noext, noext/ifact0th, 
  noext/ifact1st, pollard, pollp1, power, pp100000, squfof, 

Now we can test them all. A preliminary version of the following test indicated to me that numbers 14, 15, 16, 24, 25, and 29 were difficult to test: They either produced unrecoverable errors on my test case, or they went on way too long. So, here's the test code:

n:= nextprime(10^24)*nextprime(10^25):
for k,M in Meths do
    if k in {14,15,16,24,25,29} then next fi;
        print(k, M);
        print(CodeTools:-Usage(ifactor(n, M)))
    end try

From doing this, I couldn't find any evidence of multiple-core usage. That doesn't mean that it's not there, only that I didn't see it. Anyway, the above should give you a framework for further exploration.

A brief tutorial on "display-form" inertness in Maple

Maple offers several forms of inertness, i.e., ways to control when, if ever, expressions or particular subexpressions evaluate. I'm not aware of names for these different forms of inertness. So, for now at least, I'm going to call the topic of this thread and this article "display-form" inertness because it's primarily (although not exclusively) used to manipuate the prettyprinted displayed form of expressions while keeping them within manifestly mathematical data structures that--unlike for example strings--are only slight variations of their noninert (or active) counterparts. This form of inertness is mostly controlled by using the command value, by prefixing symbols and operators with %, and by the package InertForm. The package Typesetting and procedures whose names begin `print/ are also involved, although this aspect is usually transparent to the end user and often also transparent to the programmer who's one step removed from the end user.

The option inert=true vs. inert=false in package InertForm: The inert=true and inert=false options to some commands in InertForm do not change whether the expression is inert; they only change whether the inert operators are displayed as inert (i.e., in gray) rather than as normal (i.e., in blue).[*1] And it's only the inert operators, not the entire expression, whose color changes.

So now that I've pointed out the subtlety of this difference, do you notice it? It's clearly visible in the above Answers by @acer and @Kitonum .

If neither of these options is used, that's equivalent to inert=true.

Making symbols inert with %: Any Maple symbol (e.g., command name (if it's a procedure), function name, variable, named constant (e.g., Piinfinity)) can be made inert by prefixing an %. Although these inert forms usually "do nothing" (i.e., they are truly "inert" as in the common usage of that word), that is not actually a requirement. Executable code in procedures can be supplied for them, and those procedures can return an inert form as an unevaluated function. One possible use for that is to do syntax checking of arguments without doing the computation. An example of this is shown in the Reply below.

Special-character quoting (e.g., `%Pi`) is usually not needed (in 1D input (I don't know about 2D Input)) when prefixing %, but it is allowed. Multiple levels of inertness can be applied by prefixing multiple %. In that case, the quoting is required.

Inert arithmetic infix operators: The five basic[*2] infix arithmetic operators as well as noncommutative multiplication are available in inert forms for infix use: %^%.%*%/%+%-. These can be used infix in 1D input. I don't know about 2D Input, but your experience mentioned above suggests that they can't be. As Kitonum showed, they can always be used in prefix (aka functional) form, in which case they're just symbols, and the previous paragraph applies. You cannot create multiple levels of inertness for the infix forms.

You cannot use quotes with any infix operators (regardless of whether they're inert) if they're being used in infix form. The quotes ` ` (sometimes called "back quotes") turn anything into a symbol. Infix operators cannot be symbols, although almost all of them have an equivalent symbol form that allows them to be used in functional form (just like Kitonum used `%*` above).

Like any "multiple-punctuation-mark" infix operator in almost any computer language[*3], you cannot separate the characters with any white space. And of course ab or `ab` (those are identical) and `a b` are distinct symbols.

Other historically inert commands: The commands EvalDiffIntIntatLimitSum, Product, Normal, and Hypergeom (and perhaps I've missed some) are inert forms of the corresponding uncapitalized commands. These are much older than is using for inertness; however, their form of inertness is essentially the same as with %.

There are many other cases of two Maple commands whose names only differ by capitalization and the capitalized form is inert; yet the lowercase form is not an active form for it. In these cases, the capitalized form is usually a data container for input in some specialized domain of computation such as modmodp1, and evala. See their help pages for examples of these inert functions.

Removing inertness with value: The value command applied to an expression removes one level of inertness from every inert symbol (including the historical ones from the previous paragraph) and inert infix operator in the expression. If special techniques need to be used to remove the inertness, value has a facility for supplying procedures[*4] for that. Likewise, any otherwise unused symbol, regardless of any special characters, can be made into an inert form of any other symbol (i.e., there need not be any connection between the spellings of the inert and active forms) by teaching[*4] value about it. The previous sentence applies to symbols only, not to infix operators. It is not necessary that any noninert form produced by value actually does anything, i.e., it could also be "inert" for all practical purposes.

The InertForm package: In summary, there are numeous ways to use inertness without using the InertForm package. However, one way (among several) that the package is quite useful is for working with expressions that are so volatile that they must be entered as strings before being Parsed by the package into mathematical expressions. Examples of such volatile expressions are those that would otherwise be changed by Maple's numerous forms of automatic simplification.

The command Typeset produces displayable expressions whose detailed internal structure can be examined with lprint. This is useful for learning how the Typesetting package is used (essentially as a mark-up language).

`print/...procedures: A procedure named `print/...where ... is a function name (usually an essentially inert function, although its inertness need not be formally imposed by or by any other means) is one of the easier ways to create customized prettyprinted displays for inert functions. For example, the inert function can be broken down into simpler inert functions whose prettyprinted display is already handled by the system. An example of this is in the Reply below. 

These procedures could manipulate the output of InertForm:-Typeset. An example of this---performing a small "surgery" to correct the blue exclamation points---is also in the Reply below.

Relevant help pages: The most-important help page for more details on what I've presented here is ?value. Also relevant are ?operator,precedence, ?print, and of course ?InertForm (in particular, ?InertForm,Display and ?InertForm,Typeset). Some more details about using "punctuation mark" operators as function symbols can be found at ?use.

The commands FromInert and ToInert are not related to the "display-form" inertness being discussed in this thread and this article. They are oriented towards debugging, automatic code generation, and what I call "code surgery".

I've used several technical terms in this article without providing definitions. These appear in upright bold if they're available verbatim as Maple commands, keywords, etc., or in italics if they're being used in a technical sense (such as specific to formal computer science) rather than as common English. Feel free to ask me for specific definitions of any of these.

[*1] These styles, such as the colors gray and blue, could be changed if you really wanted. I've never seen it done, and I'd recommend against it unless you were a master craftsman at typography and document layout, lest you significantly reduce the readability of the output, especially its readability to experienced Maple users who are expecting certain styles.

An example in the Reply below corrects a bug that causes the color of the exclamation points to be default blue rather than inert gray in inert factorials.

[*2] Personally, I'd also consider mod to be a "basic" infix arithmetic operator, but it doesn't have an infix inert form. Just like most of the infix "punctuation-mark" operators, its inert infix effect can be achieved by quoting, i.e., by using `%mod` as a function symbol.

[*3] Fortran, at least in some of its older forms, is unlike most languages in that it has some exceptions to these white space restrictions.

[*4] Specifically, this is done by defining a procedure named `value/...where the ... is the inert function symbol. The details are given at ?value. An example is shown in the Reply immediately below.

There's an option to dsolve whose sole purpose is to get around this specific error. It's called continuation. See if you can find the help for it. It's probably on the help page ?dsolve,bvp_numeric,advanced (that's off the top of my head, so, not sure about that page name). If you haven't figured it out in an hour or so, I'll post an Answer. 

Also, I've answered this exact question on MaplePrimes many times. If you're lucky, you may be able to find such an Answer. Unfortunately, I can't, because the Search Tool here is horrible. It would likely take several days to hand scan through all my Answers looking for the ones related to numeric solution of boundary-layer BVPs.

Good Question; Vote Up.

There are a few snags in your formulation. These are not quite mathematical issues but rather issues of the level of symbolic abstraction, symbolic simplication, and when (if ever) variables will be given numeric values. My solution below assumes that n be specified up front (as you have done), and i be specified before the merger of the piecewises. The result for each i is a single piecewise function that is branched wrt x/h not x. So, it doesn't quite make sense that phi[i] be specified as a function of x, because this process requires x and h to be names while phi is operating. Note that the result can be expressed as a function of x/h, which for a piecewise makes it much simpler than a function of x and considered independently. The result of phi[i](x) can be applied to numeric values of x and h. I chose to do it this way so that the final result would be most amenable to further symbolic processing, such as integration. If instead the goal was a crude numeric application such as plotting, these subtler details wouldn't matter much. Let me know if you'd like any of these details rearranged.


I can't post the worksheet inline, but you can download it. Here is the code portion:

S:= x-> piecewise(
     x <= -2, 0,
     x <= -1, (2+x)^3/4,
     x <= 0,  (2+x)^3/4 - (1+x)^3,
     x <= 1,  (2-x)^3/4 - (1-x)^3,
     x <= 2,  (2-x)^3/4,
n:= 9:
_phi:= proc(x::name)
    i, xi, h, R,
    T:= table([
        0=    S(xi) - 4*S(xi+1),
        1=    S(xi-1) - S(xi+1),
        _n=   S(xi-_n) - S(xi-(_n+2)),
        _n+1= S(xi-(_n+1)) - 4*S(xi-(_n+2))
    if not procname::indexed or nops(procname) <> 1 then
        error "phi must be called with a single index"
    i:= op(procname);
    if not i::integer[0.._n+1] then
            "phi's index must be an integer from 0 to %1 but"
            " received %2", _n+1, i
    R:= `if`(assigned(T[i]), T[i], S(xi-i));    
    unapply(subs(xi= x/h, convert(R, piecewise, xi)), [x,h])
end proc
phi:= subs(_n= n, eval(_phi))
for k in [0,1,iquo(n,2),n,n+1] do print(i=k, phi[k](x)(x,h)) od:




3 4 5 6 7 8 9 Last Page 5 of 332