It's not quite clear, to me, whether you intend on using plots[odeplot] or not. I am going to suppose that you do not, and that you do intend on generating the plot data by using your Taylor4 routine.

You appear to know the value 'a' before you generate the plot data. So, pass that value of 'a' as additional data to the routine that you write which produces the plot data. At each step, after the current h is determined, have that routine examine the current t[i] and t[i+1], to see whether 'a' lies between them. When it does, have the routine generate X[t[i+1]] from t[i+1]=t[i]+h, just as usual. But also have it compute X[a] and save that as a special separate value. It could compute X[a] using a special instance of the stepsize h_i_a=a-t[i] . That special step would extend the solution from X[t[i]] to X[a]. But the rest of the whole solution up to X[t[n]] would get computed from X[t[i+1]] as normal. Then, when finished, have the routine return both the plot data (a listlist L for pointplot, or whatever) as well as X[a] as the second part of the return expression sequence.

In other words, when Taylor4 gets to t[i] the value just below 'a', have it compute both X[a] and X[t[i+1]] using two different step sizes. Save X[a]. But compute further using X[t[i+1]] as usual.

Eg, for multiple special values,

Taylor4({deq,ic},special=[a,b,c]);

could return,

L, [X[a],X[b],X[c]]

where a,b, and c are values you want treated specially in that manner.

Does that make sense? Note that you could often have your solver just continue from X[a], as if 'a' were the next natural step. Whether to do so, ignoring the usual h at that  point, depends on how you want the solver compute h. If you have a scheme that estimates h[i+1] based in some way upon h[i], then you may not want to mess up that sequence of stepsizes. a-t[i]=h_i_a may sometimes be very small, and you may well not wish to continue with such a small step size. You may wish to revert to your usual step size. (h could even be constant, as an argument to Taylor4, but I don't know whether that is so.)

There are other variants on these ideas. If 'a' is very close to t[i] then it may not be good to genarate a very small h_i_a=a-t[i] for computing X[a] = X[t[i]+h_i_a] as doing so might result in too much numerical error. It might be better to go back to t[i-1] and X[t[i-1]] and generate a special step size a-t[i-1]. That special step size would be very close to the original step size t[i]-t[i-1]. And, once again, whether to proceed onwards from X[a] or from X[t[i]] as the start of the remaining iterations would depend things like whether you'd wish to have the introduction of 'a' mess with your usual behaviour of the solver.

acer

It's not quite clear, to me, whether you intend on using plots[odeplot] or not. I am going to suppose that you do not, and that you do intend on generating the plot data by using your Taylor4 routine.

You appear to know the value 'a' before you generate the plot data. So, pass that value of 'a' as additional data to the routine that you write which produces the plot data. At each step, after the current h is determined, have that routine examine the current t[i] and t[i+1], to see whether 'a' lies between them. When it does, have the routine generate X[t[i+1]] from t[i+1]=t[i]+h, just as usual. But also have it compute X[a] and save that as a special separate value. It could compute X[a] using a special instance of the stepsize h_i_a=a-t[i] . That special step would extend the solution from X[t[i]] to X[a]. But the rest of the whole solution up to X[t[n]] would get computed from X[t[i+1]] as normal. Then, when finished, have the routine return both the plot data (a listlist L for pointplot, or whatever) as well as X[a] as the second part of the return expression sequence.

In other words, when Taylor4 gets to t[i] the value just below 'a', have it compute both X[a] and X[t[i+1]] using two different step sizes. Save X[a]. But compute further using X[t[i+1]] as usual.

Eg, for multiple special values,

Taylor4({deq,ic},special=[a,b,c]);

could return,

L, [X[a],X[b],X[c]]

where a,b, and c are values you want treated specially in that manner.

Does that make sense? Note that you could often have your solver just continue from X[a], as if 'a' were the next natural step. Whether to do so, ignoring the usual h at that  point, depends on how you want the solver compute h. If you have a scheme that estimates h[i+1] based in some way upon h[i], then you may not want to mess up that sequence of stepsizes. a-t[i]=h_i_a may sometimes be very small, and you may well not wish to continue with such a small step size. You may wish to revert to your usual step size. (h could even be constant, as an argument to Taylor4, but I don't know whether that is so.)

There are other variants on these ideas. If 'a' is very close to t[i] then it may not be good to genarate a very small h_i_a=a-t[i] for computing X[a] = X[t[i]+h_i_a] as doing so might result in too much numerical error. It might be better to go back to t[i-1] and X[t[i-1]] and generate a special step size a-t[i-1]. That special step size would be very close to the original step size t[i]-t[i-1]. And, once again, whether to proceed onwards from X[a] or from X[t[i]] as the start of the remaining iterations would depend things like whether you'd wish to have the introduction of 'a' mess with your usual behaviour of the solver.

acer

Congratulations, Alejandro. And thanks also to Axel and Alec, for their input this past month. I feel that we benefit not only from their expertise and advice but also from their constructive opinions and criticism.

acer

The round-bracket indexing for Matrix and Vectors (well, rtables) is new in Maple 12. See the ?updates,Maple12,programming help-page for some explanation of it.

It should also be in the ?updates,Maple12,language help-page, but isn't. It should also be in the ?LinearAlgebra,General,MVassignment and ?LinearAlgbera,General,MVextract help-pages, but isn't.

This illustrates why I voted for better documentation in the current poll.

acer

Ok. I guess that the OP was slightly unaware of Maple's syntax (and/or ramification) for 2-argument arctan..

Your post below, completing Alejandro's suggestion to obtain the result and apply it, is nice.

acer

Ok. I guess that the OP was slightly unaware of Maple's syntax (and/or ramification) for 2-argument arctan..

Your post below, completing Alejandro's suggestion to obtain the result and apply it, is nice.

acer

It still doesn't hold for all real A, B, and t.

> eq1 := A*sin(t) + B*cos(t) = C*sin(t+theta);
                 eq1 := A sin(t) + B cos(t) = C sin(t + theta)
 
> eq2 := theta = arctan(B/A);
                          eq2 := theta = arctan(B/A)
 
> eq3 := C = sqrt(A^2+B^2);
                                         2    2 1/2
                            eq3 := C = (A  + B )
 
> eval(eq1,[eq2,eq3]);
                                    2    2 1/2
            A sin(t) + B cos(t) = (A  + B )    sin(t + arctan(B/A))
 
> eval(eval(eq1,[eq2,eq3]),[A=-1,B=1,t=Pi]);
                                    -1 = 1

acer

It still doesn't hold for all real A, B, and t.

> eq1 := A*sin(t) + B*cos(t) = C*sin(t+theta);
                 eq1 := A sin(t) + B cos(t) = C sin(t + theta)
 
> eq2 := theta = arctan(B/A);
                          eq2 := theta = arctan(B/A)
 
> eq3 := C = sqrt(A^2+B^2);
                                         2    2 1/2
                            eq3 := C = (A  + B )
 
> eval(eq1,[eq2,eq3]);
                                    2    2 1/2
            A sin(t) + B cos(t) = (A  + B )    sin(t + arctan(B/A))
 
> eval(eval(eq1,[eq2,eq3]),[A=-1,B=1,t=Pi]);
                                    -1 = 1

acer

Are you suggesting a conversion like the following,

A*cos(t) + B*sin(t) = C*sin(t+theta)

where,

theta = arctan(B/A)
C = sqrt(A^2+B^2)

When A=1, B=-1, and t=Pi those are not equivalent, are they? Isn't it a problem, if A>0 and B<0?

acer

Are you suggesting a conversion like the following,

A*cos(t) + B*sin(t) = C*sin(t+theta)

where,

theta = arctan(B/A)
C = sqrt(A^2+B^2)

When A=1, B=-1, and t=Pi those are not equivalent, are they? Isn't it a problem, if A>0 and B<0?

acer

Jacques, do you think that the zero-tester in use by eval should be configurable?

acer

I would suggest not doing that, if avoidable. The masking of the actual entries by dummy names may lead to use of a "hidden zero" as a pivot during computation of the inverse. If that happened (though perhaps quite unlikely) then the result would be incorrect. Zeros could appear in denominators upon resubstitution and unmasking of the actual entries, with numeric exception error ensuing.

But now we know that the actual entries are not just simple names, even if they recur with a given structure. So is the huge memory use due to printing of results, or perhaps also due to computational grind during the inversion? Instead of doing the inversion using a context-menu, as was described, do it by issuing the relevant command. And suppress printing of the output, terminating that command with a colon instead of a semicolon (as alec suggested).

The code in question could also be uploaded to this site, using its File Manager.

acer

I would suggest not doing that, if avoidable. The masking of the actual entries by dummy names may lead to use of a "hidden zero" as a pivot during computation of the inverse. If that happened (though perhaps quite unlikely) then the result would be incorrect. Zeros could appear in denominators upon resubstitution and unmasking of the actual entries, with numeric exception error ensuing.

But now we know that the actual entries are not just simple names, even if they recur with a given structure. So is the huge memory use due to printing of results, or perhaps also due to computational grind during the inversion? Instead of doing the inversion using a context-menu, as was described, do it by issuing the relevant command. And suppress printing of the output, terminating that command with a colon instead of a semicolon (as alec suggested).

The code in question could also be uploaded to this site, using its File Manager.

acer

Under the Classic or commandline interfaces, one should be able to issue,

plotsetup(maplet);

Doing that makes a big difference for the output of, say, plots:-polarplot(1). Or for plots with typeset 2D Math in the caption.

A disadvantage is that control doesn't appear to return to the Classic worksheet until the pop-up plot window is dispelled.

acer

I did the above in the commandline (TTY) interface. Besides not having to do 2D Math typesetting of the results, that interface also has default support for common subexpression representation in the output.

Then I tried LinearAlgebra:-MatrixInverse(p(10,t)) in the Maple 12 Standard Java GUI. Maple can compute that result in a few seconds. But the Standard GUI is still chugging away on it after even a few minutes. And the GUI itself is using a lot of memory (500MB at this point).

Now, I don't know whether the huge amount of time and memory that the Standard GUI is using to display the result for size n=10 is due to lack of common subexpression labelling or to 2D Math output typesetting. But the TTY interface takes just a few pages to display it, with common subexpression labelling.

The Standard GUI also issued a slew of Java error messages to my console, starting with, `Exception in thread "Timer-1" java.lang.OutOfMemoryError: Java heap space'. As I type this, that example is still stalled on printing, and the Standard GUI is white (blank) and hence unresponsive. I shall have to kill the session (without being able to save the worksheet).

I tried it in Maple's Classic graphical interface. It worked beautifully.

acer