epostma

1494 Reputation

19 Badges

17 years, 59 days
Maplesoft

Social Networks and Content at Maplesoft.com

I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

MaplePrimes Activity


These are replies submitted by epostma

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.

By the way, to obtain the plot that you draw, another option would be to do this:

points := PointPlot(Y, xcoords = convert(X, list));
curve := plot(Q, x = 0 .. 115);
plots[display]([points, curve], title="Mexican Population Growth");

Of course you can set thickness, style etc. either individually for points or curves, or for both of them together by specifying the argument with the display command.

Erik.

By the way, to obtain the plot that you draw, another option would be to do this:

points := PointPlot(Y, xcoords = convert(X, list));
curve := plot(Q, x = 0 .. 115);
plots[display]([points, curve], title="Mexican Population Growth");

Of course you can set thickness, style etc. either individually for points or curves, or for both of them together by specifying the argument with the display command.

Erik.

I agree with Alec -- this is the correct solution. Thanks Robert!

Erik Postma
Maplesoft

I agree with Alec -- this is the correct solution. Thanks Robert!

Erik Postma
Maplesoft

Hi protrader,

I can see a nice way to set the number of components to a fixed value, and a different way to bound the step size to a maximum, but I'm not sure if there's an elegant way to do both. I won't post any ready-made procedures, because there's too many options currently.

If you want to fix the number of components, then don't add the numbers in between in order: just pick a number that you haven't picked before in each iteration of a loop. How to go about this depends on the required efficiency and the problem size you expect to be typical.

If you want to bound the step size, do pick the numbers in order, like Doug does. Every time pick a number in the range [1 .. your upper bound] and add that to the end of your list, unless it would be at least your upper bound minus one, in which case you are done.

Hmmm... I guess you could combine these ideas by always splitting a long interval with a random point in its interior. However, with a few unlucky picks of the random generator, you might not get a solution whereas there does exist a solution; e.g., if you want to split the interval [0..100] into four parts (so with three numbers in between), each of which is at most 26 units long, and the random number generator picks the number 99 as your first boundary, then you're toast.

Another approach would be to "seed" the process by initially picking maximum-size intervals, then later splitting them further if you have extra intervals "left over". But this would give a pretty nonrandom result for some inputs.

 

I'll think about this a bit more in the next couple of days, and if I come up with something I'll post it here.

Erik Postma
Maplesoft.

First 18 19 20 21 Page 20 of 21