dom, DD, and sp4 are defined in the loop (technically, they're probably tables, unless you've assigned something to either variable outside the loop you've shown above), so [dom[i], DD[i], sp4[i]] should be a well defined point for i in 1 .. NII. I think what you're proposing should work. Does it give an error?

Erik Postma
Maplesoft.

dom, DD, and sp4 are defined in the loop (technically, they're probably tables, unless you've assigned something to either variable outside the loop you've shown above), so [dom[i], DD[i], sp4[i]] should be a well defined point for i in 1 .. NII. I think what you're proposing should work. Does it give an error?

Erik Postma
Maplesoft.

For your first question: I don't understand why it doesn't work. If hvalues is generated by mapping over dTvalues, then they should have the same rtable dimensions. What does rtable_dims(dTvalues) return? And how about rtable_dims(hvalues)?

For your second question: do I understand correctly that you're confident that the loop iterates over the values for T and d that you want, and that you just want to generate a point plot of the generated points? I think it should work if you just type

plots:-pointplot3d(points);

after the fragment that you posted earlier.

Hope this helps, and all the best to you, too,

Erik Postma
Maplesoft.

For your first question: I don't understand why it doesn't work. If hvalues is generated by mapping over dTvalues, then they should have the same rtable dimensions. What does rtable_dims(dTvalues) return? And how about rtable_dims(hvalues)?

For your second question: do I understand correctly that you're confident that the loop iterates over the values for T and d that you want, and that you just want to generate a point plot of the generated points? I think it should work if you just type

plots:-pointplot3d(points);

after the fragment that you posted earlier.

Hope this helps, and all the best to you, too,

Erik Postma
Maplesoft.

OK - I see that dom and DD contain values in a sort of time * d space, and you're iterating over the diagonal as time goes from 0 to 10 and d goes from -30 to 30. I don't see what you're doing with all the variants of sp.

Best,

Erik Postma
Maplesoft.

OK - I see that dom and DD contain values in a sort of time * d space, and you're iterating over the diagonal as time goes from 0 to 10 and d goes from -30 to 30. I don't see what you're doing with all the variants of sp.

Best,

Erik Postma
Maplesoft.

If you know where the singular point is, say at 0, then you could put a few extra points close to it. For example,

dvalues := Vector([seq(d, d=-50..-5, 5), -1, -1/10, 1/10, 1, seq(d, d=5..50, 5)]);

For a 3d plot you could do the same as for 2d; you just have to supply values for both T and d now. So you define g and h as you're showing above, then

dTvalues := Vector([seq(seq([d, T], d=-50..50, 10), T=1/2..5, 1/2)]);
hvalues := map(dTv -> evalf[6](eval(h, [d = dTv[1], T = dTv[2]])), dTvalues);
# plotting is now slightly more tricky:
plots:-pointplot3d([seq([op(dTvalues[i]), hvalues[i]], i=rtable_dims(dTvalues))]);

However, for a 3D plot you need many more points before you have a convincing approximation of the surface than you need in 2D to get a decent approximation of a curve.

Hope this helps,

Erik Postma
Maplesoft.

NOTE: in an earlier version of this comment, there was a typo in the specification of i for the seq in the plot command.

If you know where the singular point is, say at 0, then you could put a few extra points close to it. For example,

dvalues := Vector([seq(d, d=-50..-5, 5), -1, -1/10, 1/10, 1, seq(d, d=5..50, 5)]);

For a 3d plot you could do the same as for 2d; you just have to supply values for both T and d now. So you define g and h as you're showing above, then

dTvalues := Vector([seq(seq([d, T], d=-50..50, 10), T=1/2..5, 1/2)]);
hvalues := map(dTv -> evalf[6](eval(h, [d = dTv[1], T = dTv[2]])), dTvalues);
# plotting is now slightly more tricky:
plots:-pointplot3d([seq([op(dTvalues[i]), hvalues[i]], i=rtable_dims(dTvalues))]);

However, for a 3D plot you need many more points before you have a convincing approximation of the surface than you need in 2D to get a decent approximation of a curve.

Hope this helps,

Erik Postma
Maplesoft.

NOTE: in an earlier version of this comment, there was a typo in the specification of i for the seq in the plot command.

Arrgh, and again I made a typo. That plotting line should have been:

plots:-pointplot3d([seq([op(dTvalues[i]), hvalues[i]], i=rtable_dims(dTvalues)]);

Sorry. I'll change my earlier post.

I'm not sure I'll have time to look at the worksheet you posted earlier, what with the holidays and all...

Erik Postma
Maplesoft.

Arrgh, and again I made a typo. That plotting line should have been:

plots:-pointplot3d([seq([op(dTvalues[i]), hvalues[i]], i=rtable_dims(dTvalues)]);

Sorry. I'll change my earlier post.

I'm not sure I'll have time to look at the worksheet you posted earlier, what with the holidays and all...

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.

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.