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

Up to the initial conditions it looks OK. There you'll want to say

incs := xr(0)=0, xI(0)=0, yr(0)=0, yi(0)=0, zr(0)=-0.5, zi(0)=0;

You need equations for all of them; you can't just name the initial point. Then you'll want to say either:

sol1 := dsolve(newsys union {incs});

or

sol1 := evalf(dsolve(evalf(newsys union {incs})));

to find the solution respecting the initial conditions. Now sol1 contains the solution in terms of xr and xI, etc; for sys2 you need the solution in terms of x1 etc. As Robert suggested, you can do that as

sol2 := [ x1(t)=eval(xr(t)+I*xi(t),sol1), 
          y1(t)=eval(yr(t)+I*yi(t),sol1), 
          z1(t)=eval(zr(t)+I*zi(t),sol1)];

Then finally, input sys2 the way you've showed above. Then define

sys3 := eval(sys2, sol2);

to substitute the values from sol2. Again you'll need to decompose into real and imaginary parts because of the conjugate, then add your initial conditions, etc.

Actually, why don't you consider it as one big system? So define sys1 and sys2 as you've done before, then say

bigsys := sys1 union sys2;
functions := indets(bigsys, anyfunc(identical(t)));
redefinitions := map(f -> f = cat(op(0, f), _R)(t) + I*cat(op(0,f), _I)(t), functions);
newfunctions := indets(redefinitions, anyfunc(identical(t))) minus functions;
bigsys_separated := eval(bigsys, redefinitions);
newsys := map(evalc @ Re, bigsys_separated) union map(evalc @ Im, bigsys_separated);
incs := map(`=`, eval(newfunctions, t=0), 0); ## for example -- enter your own initial conditions here
solution := dsolve(newsys union incs); ## or: evalf(dsolve(evalf(newsys union incs)));

Solving it takes a while though - I lost patience after 10 minutes. However, if you don't particularly need a symbolic result, you could also easily run a numeric IVP solver now -- just replace the last line by

solution := dsolve(newsys union incs, numeric, functions); 
## EDIT: I later realized this was a mistake -- the last argument should be newfunctions, not functions.

and read ?dsolve,numeric to find out what you can do with the result. (Hint: it's not equations you're getting back.) This is very fast, and can generate just as nice plots, if that's what you're looking for. (If so, the next help page to read is ?plots,odeplot.)

Hope this helps,

Erik Postma
Maplesoft.

Up to the initial conditions it looks OK. There you'll want to say

incs := xr(0)=0, xI(0)=0, yr(0)=0, yi(0)=0, zr(0)=-0.5, zi(0)=0;

You need equations for all of them; you can't just name the initial point. Then you'll want to say either:

sol1 := dsolve(newsys union {incs});

or

sol1 := evalf(dsolve(evalf(newsys union {incs})));

to find the solution respecting the initial conditions. Now sol1 contains the solution in terms of xr and xI, etc; for sys2 you need the solution in terms of x1 etc. As Robert suggested, you can do that as

sol2 := [ x1(t)=eval(xr(t)+I*xi(t),sol1), 
          y1(t)=eval(yr(t)+I*yi(t),sol1), 
          z1(t)=eval(zr(t)+I*zi(t),sol1)];

Then finally, input sys2 the way you've showed above. Then define

sys3 := eval(sys2, sol2);

to substitute the values from sol2. Again you'll need to decompose into real and imaginary parts because of the conjugate, then add your initial conditions, etc.

Actually, why don't you consider it as one big system? So define sys1 and sys2 as you've done before, then say

bigsys := sys1 union sys2;
functions := indets(bigsys, anyfunc(identical(t)));
redefinitions := map(f -> f = cat(op(0, f), _R)(t) + I*cat(op(0,f), _I)(t), functions);
newfunctions := indets(redefinitions, anyfunc(identical(t))) minus functions;
bigsys_separated := eval(bigsys, redefinitions);
newsys := map(evalc @ Re, bigsys_separated) union map(evalc @ Im, bigsys_separated);
incs := map(`=`, eval(newfunctions, t=0), 0); ## for example -- enter your own initial conditions here
solution := dsolve(newsys union incs); ## or: evalf(dsolve(evalf(newsys union incs)));

Solving it takes a while though - I lost patience after 10 minutes. However, if you don't particularly need a symbolic result, you could also easily run a numeric IVP solver now -- just replace the last line by

solution := dsolve(newsys union incs, numeric, functions); 
## EDIT: I later realized this was a mistake -- the last argument should be newfunctions, not functions.

and read ?dsolve,numeric to find out what you can do with the result. (Hint: it's not equations you're getting back.) This is very fast, and can generate just as nice plots, if that's what you're looking for. (If so, the next help page to read is ?plots,odeplot.)

Hope this helps,

Erik Postma
Maplesoft.

I guess easiest is

incs := {x1_R(0)=0, x1_I(0)=0, y1_R(0)=0, y1_I(0)=0, z1_R(0)=-1/2, z1_I(0)=0,
x_R(0)=0, x_I(0)=0, y_R(0)=0, y_I(0)=0, z_R(0)=-1/2, z_I(0)=0};

I used the all-zero initial conditions because they are a lot less typing: I could just use newfunctions, the names for the split parts of the time-variant functions (split into real and complex part), evaluate them at t=0 (that's eval(newfunctions, t=0); for more info, see ?eval ), and equate each of them with 0 (that's what the map(`=`, ..., 0) is for; for more info, see ?map ).

Hope this helps,

Erik Postma
Maplesoft.

I guess easiest is

incs := {x1_R(0)=0, x1_I(0)=0, y1_R(0)=0, y1_I(0)=0, z1_R(0)=-1/2, z1_I(0)=0,
x_R(0)=0, x_I(0)=0, y_R(0)=0, y_I(0)=0, z_R(0)=-1/2, z_I(0)=0};

I used the all-zero initial conditions because they are a lot less typing: I could just use newfunctions, the names for the split parts of the time-variant functions (split into real and complex part), evaluate them at t=0 (that's eval(newfunctions, t=0); for more info, see ?eval ), and equate each of them with 0 (that's what the map(`=`, ..., 0) is for; for more info, see ?map ).

Hope this helps,

Erik Postma
Maplesoft.

You're absolutely right. My mistake was that I didn't realize the parenthesis-indexing could be used with multiple indices. And although it makes sense to refer to A(5) if there are at least 5 entries in a Matrix, by virtue of C- or Fortran-ordering being specified, it does not make too much sense to extend it. But with two indices it does, of course.

Thanks Robert!

Erik Postma
Maplesoft.

You're absolutely right. My mistake was that I didn't realize the parenthesis-indexing could be used with multiple indices. And although it makes sense to refer to A(5) if there are at least 5 entries in a Matrix, by virtue of C- or Fortran-ordering being specified, it does not make too much sense to extend it. But with two indices it does, of course.

Thanks Robert!

Erik Postma
Maplesoft.

You need the former, with =, not with :=. The = operator indicates equations (mathematical statements saying that two things are the same), the := operator indicates assignments (programming statements that set a variable (or in the case of x1(0) := 0 a remember table entry - see ?remember) equal to some value in the future). What you're trying to do here is communicate to Maple the mathematical statement that, for example, xr(0) is equal to 0. You're not changing a variable.

Hope this helps,

Erik Postma
Maplesoft.

You need the former, with =, not with :=. The = operator indicates equations (mathematical statements saying that two things are the same), the := operator indicates assignments (programming statements that set a variable (or in the case of x1(0) := 0 a remember table entry - see ?remember) equal to some value in the future). What you're trying to do here is communicate to Maple the mathematical statement that, for example, xr(0) is equal to 0. You're not changing a variable.

Hope this helps,

Erik Postma
Maplesoft.

Hi Alex,

It looks like you're putting the set newsys of equations into the set of equations you're submitting to dsolve - whereas what you'll want to do is put the elements of newsys into that set of equations, not the set itself. This is in the line defining sol1. There's a couple ways you can solve that: for example as

dsolve({incs} union newsys)

or as

dsolve({incs, op(newsys)})

By the way, I'm not sure why you use evalf there -- are you looking for exact or approximate solutions? In the first case, you'd better not use evalf, and in the second case, you'll want to use ?dsolve,numeric instead.

Hope this helps,

Erik Postma
Maplesoft.

Hi Alex,

It looks like you're putting the set newsys of equations into the set of equations you're submitting to dsolve - whereas what you'll want to do is put the elements of newsys into that set of equations, not the set itself. This is in the line defining sol1. There's a couple ways you can solve that: for example as

dsolve({incs} union newsys)

or as

dsolve({incs, op(newsys)})

By the way, I'm not sure why you use evalf there -- are you looking for exact or approximate solutions? In the first case, you'd better not use evalf, and in the second case, you'll want to use ?dsolve,numeric instead.

Hope this helps,

Erik Postma
Maplesoft.

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.

First 17 18 19 20 21 Page 19 of 21