Well, if there would be a simple way then we would have found that earlier I would hope! :)

I can see a few strategies to speed this up, but I don't know if any of them will work.

1. Try more of the method=... arguments to Int. Unlikely to lead to a speedup.
2. It might be beneficial to make just one DE system, instead of the five copies needed now. You would be able to reuse the six x1_*, y1_*, z1_* variables, but you would need five different copies of x_*, y_*, and z_*, I think. This might lead to a speedup, but don't expect miracles.
3. Best would be if you could make the integral h (in real and imaginary parts, probably) part of the DE system. I did not see an easy way to do that. The technique I mean is that if in the end you want to compute int(f(t), t=0..1) then you can just add a variable F to your system defined by F(0)=0 and F'(t) = f(t). Then F(1)=int(f(t), t=0..1). This generally leads to a substantial speedup if you can pull it off. However, I don't directly see a way to do that here easily. You could rewrite h to a nested integral and write something like F'(t) = int(g, tau=0..1-t, numeric), but I'm not sure how well that will work.

Best of luck.

Erik Postma
Maplesoft.

Well, if there would be a simple way then we would have found that earlier I would hope! :)

I can see a few strategies to speed this up, but I don't know if any of them will work.

1. Try more of the method=... arguments to Int. Unlikely to lead to a speedup.
2. It might be beneficial to make just one DE system, instead of the five copies needed now. You would be able to reuse the six x1_*, y1_*, z1_* variables, but you would need five different copies of x_*, y_*, and z_*, I think. This might lead to a speedup, but don't expect miracles.
3. Best would be if you could make the integral h (in real and imaginary parts, probably) part of the DE system. I did not see an easy way to do that. The technique I mean is that if in the end you want to compute int(f(t), t=0..1) then you can just add a variable F to your system defined by F(0)=0 and F'(t) = f(t). Then F(1)=int(f(t), t=0..1). This generally leads to a substantial speedup if you can pull it off. However, I don't directly see a way to do that here easily. You could rewrite h to a nested integral and write something like F'(t) = int(g, tau=0..1-t, numeric), but I'm not sure how well that will work.

Best of luck.

Erik Postma
Maplesoft.

Do you mean that the first computation (with d = 1/2) does not lead to a result?

For the second computation, with all the different values for d, I meant that for the situation where you have defined h with a value for digits (and possibly method). If you don't then the following returns a plot for me after a couple of minutes (gathering up all the necessary pieces):

restart;
sys1:={diff(x1(t),t)=(3+5*I)*x1(t)-sqrt(2)*(cos(3)
  +I*sin(3))*y1(t)-2*I*z1(t),diff(y1(t),t)=
  (3-5*I)*y1(t)-sqrt(2)*(cos(3)-I*sin(3))*x1(t)+2*I*z1(t),diff(z1(t),t)=0.25-0.3*z1(t)-2*I*(x1(t)-y1(t))}:
sys2:={diff(x(t),t)=(3+5*I)*x(t)-sqrt(2)*(cos(3)+I*sin(3))*y(t)-2*I*z(t)-2*I*(cos(5*t)-1)*z1(t),
  diff(y(t),t)=conjugate(x(t)),diff(z(t),t)=0.25-0.3*z(t)-2*I*(x(t)-y(t))-2*I*(cos(5*t)-1)*(x1(t)-y1(t))}:
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);
constant_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_I(0)=0, y_I(0)=0, z_I(0)=0}:
S1:= dsolve(newsys union constant_incs union {x_R(0)=1, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
S2:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=1, z_R(0)=0},
numeric, output=listprocedure):
S3:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=1},
numeric, output=listprocedure):
S4:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
f4 := eval(x_R(t) + I*x_I(t), S4);
f1 := eval(x_R(t) + I*x_I(t), S1) - f4;
f2 := eval(x_R(t) + I*x_I(t), S2) - f4;
f3 := eval(x_R(t) + I*x_I(t), S3) - f4;
S5 := dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=-1/2}, numeric, output=listprocedure):
z_R_5 := eval(z_R(t), S5);
y_5 := eval(y_R(t) + I*y_I(t), S5);
g := 2*Re(exp(-I*d*tau)*((1/2+z_R_5(t))*f1(tau)+f4(tau)-f3(tau)/2*y_5(t)));
h := Int(g, [tau = 0..1-t, t=0..1]);
dvalues := Vector([seq(d/10, d=0..20)]);
hvalues := map(dv -> evalf[6](eval(h, d = dv)), dvalues);

Now at this point, it turns out that for d = 1.3, Maple cannot evaluate the integral with the requested precision for some reason, so the integral returns unevaluated. But if we drop one digit it does work and leads to a plot:

hvalues[14] := evalf[5](hvalues[14]);
plots:-pointplot(dvalues, hvalues);

Hope this helps,

Erik Postma
Maplesoft.

Do you mean that the first computation (with d = 1/2) does not lead to a result?

For the second computation, with all the different values for d, I meant that for the situation where you have defined h with a value for digits (and possibly method). If you don't then the following returns a plot for me after a couple of minutes (gathering up all the necessary pieces):

restart;
sys1:={diff(x1(t),t)=(3+5*I)*x1(t)-sqrt(2)*(cos(3)
  +I*sin(3))*y1(t)-2*I*z1(t),diff(y1(t),t)=
  (3-5*I)*y1(t)-sqrt(2)*(cos(3)-I*sin(3))*x1(t)+2*I*z1(t),diff(z1(t),t)=0.25-0.3*z1(t)-2*I*(x1(t)-y1(t))}:
sys2:={diff(x(t),t)=(3+5*I)*x(t)-sqrt(2)*(cos(3)+I*sin(3))*y(t)-2*I*z(t)-2*I*(cos(5*t)-1)*z1(t),
  diff(y(t),t)=conjugate(x(t)),diff(z(t),t)=0.25-0.3*z(t)-2*I*(x(t)-y(t))-2*I*(cos(5*t)-1)*(x1(t)-y1(t))}:
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);
constant_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_I(0)=0, y_I(0)=0, z_I(0)=0}:
S1:= dsolve(newsys union constant_incs union {x_R(0)=1, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
S2:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=1, z_R(0)=0},
numeric, output=listprocedure):
S3:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=1},
numeric, output=listprocedure):
S4:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
f4 := eval(x_R(t) + I*x_I(t), S4);
f1 := eval(x_R(t) + I*x_I(t), S1) - f4;
f2 := eval(x_R(t) + I*x_I(t), S2) - f4;
f3 := eval(x_R(t) + I*x_I(t), S3) - f4;
S5 := dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=-1/2}, numeric, output=listprocedure):
z_R_5 := eval(z_R(t), S5);
y_5 := eval(y_R(t) + I*y_I(t), S5);
g := 2*Re(exp(-I*d*tau)*((1/2+z_R_5(t))*f1(tau)+f4(tau)-f3(tau)/2*y_5(t)));
h := Int(g, [tau = 0..1-t, t=0..1]);
dvalues := Vector([seq(d/10, d=0..20)]);
hvalues := map(dv -> evalf[6](eval(h, d = dv)), dvalues);

Now at this point, it turns out that for d = 1.3, Maple cannot evaluate the integral with the requested precision for some reason, so the integral returns unevaluated. But if we drop one digit it does work and leads to a plot:

hvalues[14] := evalf[5](hvalues[14]);
plots:-pointplot(dvalues, hvalues);

Hope this helps,

Erik Postma
Maplesoft.

What do you mean by "didn't get the value"? h is supposed to be an unevaluated integral; that's necessary because there's still a variable d in there. This is why you need evalf(eval(h, d = 1/2)) or some other value in order to obtain a number.

If this did for some reason not work, please let us know what the error message or result was. Also, have you tried this approach without the method=... argument, so that you can use the original ranges and don't need the piecewise?

Hope this helps,

Erik Postma
Maplesoft.

What do you mean by "didn't get the value"? h is supposed to be an unevaluated integral; that's necessary because there's still a variable d in there. This is why you need evalf(eval(h, d = 1/2)) or some other value in order to obtain a number.

If this did for some reason not work, please let us know what the error message or result was. Also, have you tried this approach without the method=... argument, so that you can use the original ranges and don't need the piecewise?

Hope this helps,

Erik Postma
Maplesoft.

method=_MonteCarlo would be one of the options, yes, but it will only work for low accuracy (as the help page ?evalf,Int says) and additionally it will only work for a rectangular region. You could transform the problem to something like

h := Int(piecewise(tau > 1-t, 0, g), [tau=0..1, t=0..1], digits=5, method=_MonteCarlo));
evalf(eval(h, d=1/2));

but I don't think it will necessarily be the most efficient. The _cuhre method has that same issue (needs a rectangular region of integration), but it might be useful to try the same trick there, especially if you need higher accuracy. Otherwise Maple will use nested 1D solving; you can experiment with the other methods to see if there's one that happens to be fast than the default for your problem, but maybe the standard heuristics are best in that case.

In general I would agree that a good way to obtain a plot for this problem, where computation of points is so expensive, is to explicitly compute a number of points; once you've found a definition of h that you're happy with (which should not involve evalf, since there's still a variable d in there), you could do something like this:

dvalues := Vector([seq(d/10, d=0..20)]);
hvalues := map(dv -> evalf(eval(h, d = dv)), dvalues);
plots:-pointplot(dvalues, hvalues);

On the other hand, if you need to do the plot just once, it might be worth it to wait a few more minutes and instead define

hfunction := dv -> evalf(eval(h, d = dv));
plot(hfunction, 0..2);

or whatever your valid range is.

Hope this helps,

Erik Postma
Maplesoft.

method=_MonteCarlo would be one of the options, yes, but it will only work for low accuracy (as the help page ?evalf,Int says) and additionally it will only work for a rectangular region. You could transform the problem to something like

h := Int(piecewise(tau > 1-t, 0, g), [tau=0..1, t=0..1], digits=5, method=_MonteCarlo));
evalf(eval(h, d=1/2));

but I don't think it will necessarily be the most efficient. The _cuhre method has that same issue (needs a rectangular region of integration), but it might be useful to try the same trick there, especially if you need higher accuracy. Otherwise Maple will use nested 1D solving; you can experiment with the other methods to see if there's one that happens to be fast than the default for your problem, but maybe the standard heuristics are best in that case.

In general I would agree that a good way to obtain a plot for this problem, where computation of points is so expensive, is to explicitly compute a number of points; once you've found a definition of h that you're happy with (which should not involve evalf, since there's still a variable d in there), you could do something like this:

dvalues := Vector([seq(d/10, d=0..20)]);
hvalues := map(dv -> evalf(eval(h, d = dv)), dvalues);
plots:-pointplot(dvalues, hvalues);

On the other hand, if you need to do the plot just once, it might be worth it to wait a few more minutes and instead define

hfunction := dv -> evalf(eval(h, d = dv));
plot(hfunction, 0..2);

or whatever your valid range is.

Hope this helps,

Erik Postma
Maplesoft.

This is what I did:

S5 := dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=-1/2}, numeric, output=listprocedure):
z_R_5 := eval(z_R(t), S5);
y_5 := eval(y_R(t) + I*y_I(t), S5);
g := 2*Re(exp(-I*d*tau)*((1/2+z_R_5(t))*f1(tau)+f4(tau)-f3(tau)/2*y_5(t)));
h := Int(g, [tau = 0..1-t, t=0..1]);
evalf[6](eval(h, d=1/2));

The 6 in that last command indicates that you request 6 digits of accuracy; the answer is 0.343701. (As you see, this is for d=1/2.) This took ~30 seconds on my old, slow laptop.

If your goal is to get a plot of the result as a function of d, then you'll need quite a bit of computation time. I'd strongly recommend playing with the method argument as described in ?evalf,int to see what works best for this particular problem.

Hope this helps,

Erik Postma
Maplesoft.
 

This is what I did:

S5 := dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=-1/2}, numeric, output=listprocedure):
z_R_5 := eval(z_R(t), S5);
y_5 := eval(y_R(t) + I*y_I(t), S5);
g := 2*Re(exp(-I*d*tau)*((1/2+z_R_5(t))*f1(tau)+f4(tau)-f3(tau)/2*y_5(t)));
h := Int(g, [tau = 0..1-t, t=0..1]);
evalf[6](eval(h, d=1/2));

The 6 in that last command indicates that you request 6 digits of accuracy; the answer is 0.343701. (As you see, this is for d=1/2.) This took ~30 seconds on my old, slow laptop.

If your goal is to get a plot of the result as a function of d, then you'll need quite a bit of computation time. I'd strongly recommend playing with the method argument as described in ?evalf,int to see what works best for this particular problem.

Hope this helps,

Erik Postma
Maplesoft.
 

As you can see in the example that I gave in my previous comment, f1 etc. are not expressions but procedures. That means that you can't substitute a value into it - you just apply it to a value. So you would define g sort of like this:

g:=(d,T)->2/T*(Re(exp(-I*d*tau)*((0.5+z_R(t))*f1(tau)+f4(tau)-0.5*f3(tau)*y(t)))):

Of course, you'd also have to substitute the correct procedures for z_R and y. Furthermore, I'm not sure why g is a function of d and T but not of tau, but I trust you'll have your reasons.

Hope this helps,

Erik Postma
Maplesoft.

As you can see in the example that I gave in my previous comment, f1 etc. are not expressions but procedures. That means that you can't substitute a value into it - you just apply it to a value. So you would define g sort of like this:

g:=(d,T)->2/T*(Re(exp(-I*d*tau)*((0.5+z_R(t))*f1(tau)+f4(tau)-0.5*f3(tau)*y(t)))):

Of course, you'd also have to substitute the correct procedures for z_R and y. Furthermore, I'm not sure why g is a function of d and T but not of tau, but I trust you'll have your reasons.

Hope this helps,

Erik Postma
Maplesoft.

When you claim that x(t) = x(0)*f1(t) + y(0)*f2(t) + z(0)*f3(t) + f4(t), then I guess you assume that x1(0), y1(0), z1(0) are constants - in particular that they are 0, 0, -1/2 as before? Then I think what you'll want to do is this:

constant_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_I(0)=0, y_I(0)=0, z_I(0)=0}:
S1:= dsolve(newsys union constant_incs union {x_R(0)=1, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
S2:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=1, z_R(0)=0},
numeric, output=listprocedure):
S3:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=1},
numeric, output=listprocedure):
S4:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
f4 := eval(x_R(t) + I*x_I(t), S4);
f1 := eval(x_R(t) + I*x_I(t), S1) - f4;
f2 := eval(x_R(t) + I*x_I(t), S2) - f4;
f3 := eval(x_R(t) + I*x_I(t), S3) - f4;

Now, assuming that indeed the general solution can be written as x(t) = x(0)*f1(t) + y(0)*f2(t) + z(0)*f3(t) + f4(t), we have correctly identified f1 and f4. You can now, for example, type f1(1) to obtain the answer  9.41518856579549 - 9.82217797194972 I.

Hope this helps,

Erik Postma
Maplesoft.

PS: I briefly had a different version of this response online which contained a mistake (it used x1_R(t) + I*x1_I(t) to define f1 up to f4) - if you saw it, please disregard and use this instead.

When you claim that x(t) = x(0)*f1(t) + y(0)*f2(t) + z(0)*f3(t) + f4(t), then I guess you assume that x1(0), y1(0), z1(0) are constants - in particular that they are 0, 0, -1/2 as before? Then I think what you'll want to do is this:

constant_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_I(0)=0, y_I(0)=0, z_I(0)=0}:
S1:= dsolve(newsys union constant_incs union {x_R(0)=1, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
S2:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=1, z_R(0)=0},
numeric, output=listprocedure):
S3:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=1},
numeric, output=listprocedure):
S4:= dsolve(newsys union constant_incs union {x_R(0)=0, y_R(0)=0, z_R(0)=0},
numeric, output=listprocedure):
f4 := eval(x_R(t) + I*x_I(t), S4);
f1 := eval(x_R(t) + I*x_I(t), S1) - f4;
f2 := eval(x_R(t) + I*x_I(t), S2) - f4;
f3 := eval(x_R(t) + I*x_I(t), S3) - f4;

Now, assuming that indeed the general solution can be written as x(t) = x(0)*f1(t) + y(0)*f2(t) + z(0)*f3(t) + f4(t), we have correctly identified f1 and f4. You can now, for example, type f1(1) to obtain the answer  9.41518856579549 - 9.82217797194972 I.

Hope this helps,

Erik Postma
Maplesoft.

PS: I briefly had a different version of this response online which contained a mistake (it used x1_R(t) + I*x1_I(t) to define f1 up to f4) - if you saw it, please disregard and use this instead.

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.