Hi Tom,

Interesting that the code is not working. I just checked and saw that Maple 11 and earlier do not accept vectors as arguments to min and max, so if you're running that version you might want to change that plot command to

(**) p2 := plot(result, x=min(op(convert(xd, list))) .. max(op(convert(xd, list))));

I'm not sure what you mean by "actually it isn't getting defined as a fit". As Georgios points out below, I use the long version of the commands (*), so you don't need to activate packages before using the commands. Although... if you're using a Maple version before Maple 10, you won't have the Statistics package; is that the case?

About learning Maple: I don't have any really great tips on where to find the best maple learning materials, although you could try ?portal. Myself, I learned it mostly through the manuals that come in the box, but I had the advantage of having a good number of the world's Maple experts a couple of desks away while doing it :)

Hope this helps,

Erik Postma
Maplesoft.

(*): I did it inconsistently, though, mixing table access and module access syntax. Well, either will work.

Hi Tom,

Interesting that the code is not working. I just checked and saw that Maple 11 and earlier do not accept vectors as arguments to min and max, so if you're running that version you might want to change that plot command to

(**) p2 := plot(result, x=min(op(convert(xd, list))) .. max(op(convert(xd, list))));

I'm not sure what you mean by "actually it isn't getting defined as a fit". As Georgios points out below, I use the long version of the commands (*), so you don't need to activate packages before using the commands. Although... if you're using a Maple version before Maple 10, you won't have the Statistics package; is that the case?

About learning Maple: I don't have any really great tips on where to find the best maple learning materials, although you could try ?portal. Myself, I learned it mostly through the manuals that come in the box, but I had the advantage of having a good number of the world's Maple experts a couple of desks away while doing it :)

Hope this helps,

Erik Postma
Maplesoft.

(*): I did it inconsistently, though, mixing table access and module access syntax. Well, either will work.

You can even enter it in 2D mode by just hitting the keys "x", "^", "~", "2", "right arrow" (to get out of the exponent again) and the tilde will be correctly interpreted. I'd say this is the recommended approach, when using Maple 13. (The ?elementwise operators were not available in Maple 12 or earlier.)

The code snippet pirolafire shows below doesn't make much sense to me, to be honest, and it doesn't seem to work on my machine. My guess would be that either j happened to be defined to be 1 in that session, or there was a typo in the definition of f. In the first case, it's not obvious why it works, but by digging a little we can see what's going on. The Matrix definition would expand to

Matrix(5, 1, 5, (i) -> c^2 * s[i, j]^2 - x[i, j]^2 - y[i, j]^2 - z[i, j]^2);

This means: create a five-by-one matrix, set all entries to five, then take the procedure i -> c^2*s[i,j]^2 - ... and feed it all "coordinate pairs" for the entries in the matrix, overwriting all the fives just put in there. Because of the way Maple evaluates function calls, calling that procedure with two coordinates just assigns the first coordinate (the row index) to i and discards the column coordinate. In order for this to give the correct result, j has to be defined externally to be equal to one.

In the second case, this definition of f would make more sense to me:

f := (i, j) -> c^2 * s[i, j]^2 - x[i, j]^2 - y[i, j]^2 - z[i, j]^2;

That said, I still prefer the elementwise operators. It's much cleaner and allows you to express what you want to do in a more mathematical way.

Hope this helps,

Erik Postma
Maplesoft.

You can even enter it in 2D mode by just hitting the keys "x", "^", "~", "2", "right arrow" (to get out of the exponent again) and the tilde will be correctly interpreted. I'd say this is the recommended approach, when using Maple 13. (The ?elementwise operators were not available in Maple 12 or earlier.)

The code snippet pirolafire shows below doesn't make much sense to me, to be honest, and it doesn't seem to work on my machine. My guess would be that either j happened to be defined to be 1 in that session, or there was a typo in the definition of f. In the first case, it's not obvious why it works, but by digging a little we can see what's going on. The Matrix definition would expand to

Matrix(5, 1, 5, (i) -> c^2 * s[i, j]^2 - x[i, j]^2 - y[i, j]^2 - z[i, j]^2);

This means: create a five-by-one matrix, set all entries to five, then take the procedure i -> c^2*s[i,j]^2 - ... and feed it all "coordinate pairs" for the entries in the matrix, overwriting all the fives just put in there. Because of the way Maple evaluates function calls, calling that procedure with two coordinates just assigns the first coordinate (the row index) to i and discards the column coordinate. In order for this to give the correct result, j has to be defined externally to be equal to one.

In the second case, this definition of f would make more sense to me:

f := (i, j) -> c^2 * s[i, j]^2 - x[i, j]^2 - y[i, j]^2 - z[i, j]^2;

That said, I still prefer the elementwise operators. It's much cleaner and allows you to express what you want to do in a more mathematical way.

Hope this helps,

Erik Postma
Maplesoft.

OK, so now it seems we are getting to the heart of the matter.

The differential equation you showed is about t(x), that is, a function of one variable. If you give a differential equation for a function of one variable that involves its second derivative but no higher derivatives, and if the DE is sufficiently "well behaved", then you have two degrees of freedom left. So you can give two independent equations involving the function t and its derivative at a well-defined point. In the worksheet that you posted, this point was x=-2. We typically call such equations the "initial conditions" for the DE. You could also use equations involving t and its derivatives at different points, for example t(-2) and t'(0) (in which case you would talk about boundary conditions instead of initial conditions), but you can't  give an equation involving t(x) at general values x (or at all values x). That would be infinitely many extra equations: one for each value of x. There simply is not enough "freedom" left for that after you have given your differential equation.

There are other types of differential equations where you can give boundary values at general values of an independent variable, but then you need to go to partial differential equations instead of ordinary differential equations. That means that the "subject" of the DE is a function of two or more variables; for example, f(x, y). ?pdsolve,numeric has a few examples.

So now back to your question. You have an ordinary differential equation involving the second derivative with one parameter; the one parameter gives you one additional degree of freedom, for a total of three. You now need to choose which three properties you'd like to impose on the result.

Let us know if this helps. If you're still convinced that some property should hold for all values of x, let us know why you think such a solution should exist. Then maybe we can find out together what's going on.

Hope this helps,

Erik Postma
Maplesoft.

OK, so now it seems we are getting to the heart of the matter.

The differential equation you showed is about t(x), that is, a function of one variable. If you give a differential equation for a function of one variable that involves its second derivative but no higher derivatives, and if the DE is sufficiently "well behaved", then you have two degrees of freedom left. So you can give two independent equations involving the function t and its derivative at a well-defined point. In the worksheet that you posted, this point was x=-2. We typically call such equations the "initial conditions" for the DE. You could also use equations involving t and its derivatives at different points, for example t(-2) and t'(0) (in which case you would talk about boundary conditions instead of initial conditions), but you can't  give an equation involving t(x) at general values x (or at all values x). That would be infinitely many extra equations: one for each value of x. There simply is not enough "freedom" left for that after you have given your differential equation.

There are other types of differential equations where you can give boundary values at general values of an independent variable, but then you need to go to partial differential equations instead of ordinary differential equations. That means that the "subject" of the DE is a function of two or more variables; for example, f(x, y). ?pdsolve,numeric has a few examples.

So now back to your question. You have an ordinary differential equation involving the second derivative with one parameter; the one parameter gives you one additional degree of freedom, for a total of three. You now need to choose which three properties you'd like to impose on the result.

Let us know if this helps. If you're still convinced that some property should hold for all values of x, let us know why you think such a solution should exist. Then maybe we can find out together what's going on.

Hope this helps,

Erik Postma
Maplesoft.

Whoops -- I'm usually very careful to test what I write on MaplePrimes in the most recent released version of Maple, as opposed to the version that's in development that we use for day-to-day work internally. But this time I forgot, apparently :)

Acer's approach is better anyway; I guess I just wrote zerovalue - 8 instead of (t -> zerovalue(t) - 8) out of laziness.

Erik Postma
Maplesoft.

Whoops -- I'm usually very careful to test what I write on MaplePrimes in the most recent released version of Maple, as opposed to the version that's in development that we use for day-to-day work internally. But this time I forgot, apparently :)

Acer's approach is better anyway; I guess I just wrote zerovalue - 8 instead of (t -> zerovalue(t) - 8) out of laziness.

Erik Postma
Maplesoft.

And of course I could have written that last procedure more intuitively as

zerovalue := proc(cvalue)
global solution, t, x;
  solution('parameters' = [cvalue]);
  try
    return eval(t(x), solution(0));
  catch:
    return 0;
  end try;
end proc;

Anyway. If you wanted to find the value of c where the maximum is at zero, you could use this:

zeroderiv := proc(cvalue)
global solution, t, x;
  solution('parameters' = [cvalue]);
  try
    return eval(diff(t(x), x), solution(0));
  catch:
    return -10;
  end try;
end proc;
fsolve(zeroderiv, 190 .. 210);

This returns 203.9270540, with a maximum value of about 13.

I also tried with higher settings of Digits. One thing to note is that you then need to rerun the line defining solution, otherwise it still runs with the same low Digits setting. Unfortunately, fsolve is then too clever to do this; I guess zeroderiv is numerically unstable. So I had to make do with the following loop:

Digits := 50:
solution := dsolve({ode, ics}, numeric, t(x), parameters=[c]);
high := 210.;
low := 190.;
while high - low > 10^(-40) do
  mid := low + (high - low)/2;
  if zeroderiv(mid) > 0 then
    high := mid;
  else
    low := mid;
  end if;
end do;

Hope this helps,

Erik Postma
Maplesoft

And of course I could have written that last procedure more intuitively as

zerovalue := proc(cvalue)
global solution, t, x;
  solution('parameters' = [cvalue]);
  try
    return eval(t(x), solution(0));
  catch:
    return 0;
  end try;
end proc;

Anyway. If you wanted to find the value of c where the maximum is at zero, you could use this:

zeroderiv := proc(cvalue)
global solution, t, x;
  solution('parameters' = [cvalue]);
  try
    return eval(diff(t(x), x), solution(0));
  catch:
    return -10;
  end try;
end proc;
fsolve(zeroderiv, 190 .. 210);

This returns 203.9270540, with a maximum value of about 13.

I also tried with higher settings of Digits. One thing to note is that you then need to rerun the line defining solution, otherwise it still runs with the same low Digits setting. Unfortunately, fsolve is then too clever to do this; I guess zeroderiv is numerically unstable. So I had to make do with the following loop:

Digits := 50:
solution := dsolve({ode, ics}, numeric, t(x), parameters=[c]);
high := 210.;
low := 190.;
while high - low > 10^(-40) do
  mid := low + (high - low)/2;
  if zeroderiv(mid) > 0 then
    high := mid;
  else
    low := mid;
  end if;
end do;

Hope this helps,

Erik Postma
Maplesoft

A useful but little known feature of Arrays (and Vectors, but not Matrices I think) is that you can also use them quite efficiently if you don't know the size beforehand:

(**) a := Array(1..1);
                                   a := [0]

(**) for i while rand() > 10^6 do
(**) a(i) := i:
(**) end do:
memory used=1505.0MB, alloc=662.0MB, time=28.47
memory used=1581.3MB, alloc=662.0MB, time=29.75
memory used=1657.6MB, alloc=662.0MB, time=30.98
memory used=1733.8MB, alloc=662.0MB, time=32.32
memory used=1810.1MB, alloc=662.0MB, time=33.66
memory used=1886.4MB, alloc=662.0MB, time=34.84
(**) a;
                           [ 1..1216951 1-D Array ]
                           [ Data Type: anything  ]
                           [ Storage: rectangular ]
                           [ Order: Fortran_order ]

This runs the loop a random number of times (1216951 in this case) and grows the Array whenever necessary. For us users the Array seems to grow by one element every time, but behind the scenes the kernel uses a much more efficient strategy, growing it by a certain factor if necessary. Note that this only works with  "programmer indexing", as explained on ?rtable_indexing.

Of course if you do know the number of elements you're going to need beforehand, then it's more efficient to allocate the right amount of space. Even if you know the right order of magnitude it's better to pre-allocate something of that order of magnitude. For the example above, I could have started with an Array of size 10^6, since that's the expected Array size I would end up with; but I thought starting from size 1 would illustrate the point better.

Hope this helps,

Erik Postma
Maplesoft.

A useful but little known feature of Arrays (and Vectors, but not Matrices I think) is that you can also use them quite efficiently if you don't know the size beforehand:

(**) a := Array(1..1);
                                   a := [0]

(**) for i while rand() > 10^6 do
(**) a(i) := i:
(**) end do:
memory used=1505.0MB, alloc=662.0MB, time=28.47
memory used=1581.3MB, alloc=662.0MB, time=29.75
memory used=1657.6MB, alloc=662.0MB, time=30.98
memory used=1733.8MB, alloc=662.0MB, time=32.32
memory used=1810.1MB, alloc=662.0MB, time=33.66
memory used=1886.4MB, alloc=662.0MB, time=34.84
(**) a;
                           [ 1..1216951 1-D Array ]
                           [ Data Type: anything  ]
                           [ Storage: rectangular ]
                           [ Order: Fortran_order ]

This runs the loop a random number of times (1216951 in this case) and grows the Array whenever necessary. For us users the Array seems to grow by one element every time, but behind the scenes the kernel uses a much more efficient strategy, growing it by a certain factor if necessary. Note that this only works with  "programmer indexing", as explained on ?rtable_indexing.

Of course if you do know the number of elements you're going to need beforehand, then it's more efficient to allocate the right amount of space. Even if you know the right order of magnitude it's better to pre-allocate something of that order of magnitude. For the example above, I could have started with an Array of size 10^6, since that's the expected Array size I would end up with; but I thought starting from size 1 would illustrate the point better.

Hope this helps,

Erik Postma
Maplesoft.

Thanks Schivnorr and jakubi; I'll enter both into our system.

Erik Postma
Maplesoft.

Hi ndattani,

The elementwise operator, ~, was added in Maple 13. In Maple 12, you can achieve similar results using the map command (as Doug indicated).  So instead of A~(S), you would use map(A, S). The command involving myprint2 would look something like

myprint2( Vector( map(f, S) ) );

Hope this helps,

Erik Postma
Maplesoft.

Hi ndattani,

The elementwise operator, ~, was added in Maple 13. In Maple 12, you can achieve similar results using the map command (as Doug indicated).  So instead of A~(S), you would use map(A, S). The command involving myprint2 would look something like

myprint2( Vector( map(f, S) ) );

Hope this helps,

Erik Postma
Maplesoft.