Hi Alex,

I think the easiest way to obtain f1 is to supply two sets of initial conditions: one where x(0) = y(0) = z(0) = 0 - this should give you f4(t) as the solution for x(t) - and one where x(0) = 1, y(0)=z(0)=0 - this should give you f1(t) + f4(t). Subtract the former from the latter and you're done.

Hope this helps,

Erik Postma
Maplesoft.

Hi Alex,

I can see that something like what you're posting would work if you include the option output=listprocedure in the dsolve(..., numeric) call, but using plots[odeplot] is easier, more efficient, more powerful, and potentially more accurate (in that you can make it do adaptive plotting). I'll show below.

Easier: here's one command to do what you are asking - plot x_R, z_R, and x_I (against t, presumably):

plots[odeplot](solution, [[t, x_R(t)], [t, z_R(t)], [t, x_I(t)]], 0..5);

More efficient: it doesn't use separate storage for the intermediate arrays that you're using. (Word of advice: in general, the more modern Array is to be preferred over the older array.)

More powerful: you can just as easily make a phase plot using this approach, or indeed plot any function of the variables against any other function. Even a 3D curve is possible. For example,

plots[odeplot](solution, [x_R(t), x_I(t)], 0..5, numpoints=500); # a 2D phase plot of the complex values of x
plots[odeplot](solution, [t, x_R(t), x_I(t)], 0..5, numpoints=500); # a 3D plot of the complex values of x over time
plots[odeplot](solution, [[t, arctan(x_I(t), x_R(t))], [t, ln(sqrt(x_R(t)^2 + x_I(t)^2))]], 0..5, numpoints=500); 
# a plot of the complex argument and logarithm of the norm of x over time

Potentially more accurate: I won't go into that here, but see the ?plots,odeplot help page, in particular the last example.

I'm not sure I understand your second question - why do you think the solution for x(t) can be expressed that way?

Hope this helps,

Erik Postma
Maplesoft.

Hi Alex,

I can see that something like what you're posting would work if you include the option output=listprocedure in the dsolve(..., numeric) call, but using plots[odeplot] is easier, more efficient, more powerful, and potentially more accurate (in that you can make it do adaptive plotting). I'll show below.

Easier: here's one command to do what you are asking - plot x_R, z_R, and x_I (against t, presumably):

plots[odeplot](solution, [[t, x_R(t)], [t, z_R(t)], [t, x_I(t)]], 0..5);

More efficient: it doesn't use separate storage for the intermediate arrays that you're using. (Word of advice: in general, the more modern Array is to be preferred over the older array.)

More powerful: you can just as easily make a phase plot using this approach, or indeed plot any function of the variables against any other function. Even a 3D curve is possible. For example,

plots[odeplot](solution, [x_R(t), x_I(t)], 0..5, numpoints=500); # a 2D phase plot of the complex values of x
plots[odeplot](solution, [t, x_R(t), x_I(t)], 0..5, numpoints=500); # a 3D plot of the complex values of x over time
plots[odeplot](solution, [[t, arctan(x_I(t), x_R(t))], [t, ln(sqrt(x_R(t)^2 + x_I(t)^2))]], 0..5, numpoints=500); 
# a plot of the complex argument and logarithm of the norm of x over time

Potentially more accurate: I won't go into that here, but see the ?plots,odeplot help page, in particular the last example.

I'm not sure I understand your second question - why do you think the solution for x(t) can be expressed that way?

Hope this helps,

Erik Postma
Maplesoft.

Whoops - that should be

solution := dsolve(newsys union incs, numeric, newfunctions);

Sorry! I added a note to my earlier comment, for those who revisit this later.

Erik Postma
Maplesoft.

Whoops - that should be

solution := dsolve(newsys union incs, numeric, newfunctions);

Sorry! I added a note to my earlier comment, for those who revisit this later.

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.

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.