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.

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.