Hi Mario,

I would like to experiment a little more. Would it be possible for you to upload your later revision of the code?

ps. In an earlier post in this thread, you wrote that you wanted to do something like,

p:=Taylor4(......)
p(a)
plots[odeplot](p)

It was from that that I misunderstood, and got the idea that you would know point 'a' before viewing the plot. I now see that it is the other way around: you first view the plot, and only then (due to some qualitative aspect of the plot, perhaps) know the point 'a' at which you then also want the computation p(a).

There may still be hope. Having computed all the [t[i],X[t[i]]] pairs in the data set, it should be possible to find the largest t[i] smaller than 'a' (or second largest, if the largest is too close to 'a'.) Maybe one could then use that X[t[i]] as a new "initial" point and 'a' as new final point.

acer

Hi Mario,

I would like to experiment a little more. Would it be possible for you to upload your later revision of the code?

ps. In an earlier post in this thread, you wrote that you wanted to do something like,

p:=Taylor4(......)
p(a)
plots[odeplot](p)

It was from that that I misunderstood, and got the idea that you would know point 'a' before viewing the plot. I now see that it is the other way around: you first view the plot, and only then (due to some qualitative aspect of the plot, perhaps) know the point 'a' at which you then also want the computation p(a).

There may still be hope. Having computed all the [t[i],X[t[i]]] pairs in the data set, it should be possible to find the largest t[i] smaller than 'a' (or second largest, if the largest is too close to 'a'.) Maybe one could then use that X[t[i]] as a new "initial" point and 'a' as new final point.

acer

Ok, the general idea of what I was saying should hold for your particular code. The actual implementation would have to accomodate the fact that you use a slightly different scheme, according to whether T[i]<=bd. But what I wrote above should work.

Some other advice, if I may:

• Use `add` instead of `sum` in that routine, where you add up a finite number of values. It is faster.
• Instead of all those `subs` calls, you should be able to do it faster with `eval`. Now, I remember that you asked about it elsewhere, and that you do need the successive substitutions, in order. But why not substitute into the S[i-1] (to get "new" S[i-1] formulae) outside of the 'i' loop? That is to say, get your hands on fully explicit algebraic substitution formulae, newS[i],i=1..5, just once outside the loops. And then use just a single call to eval(X[i]*F||i(t),[..newS[i]=..] inside the loop. It would run much faster, I suspect. And the routine is slow.
• Remove N as a declared local and the first line which assigns to N using n, and change the parameter n::integer to N::integer:=1000.
• I would write procedures in 1D Maple input, not 2D Math. It looks bad in 2D. And I can actually see the delay which the GUI parses and displays the code. And it is wholly unclear visually whether the subscripted names are table references or not.
• Use the parameter names instead of args[i], all over, for consistency.

acer

Ok, the general idea of what I was saying should hold for your particular code. The actual implementation would have to accomodate the fact that you use a slightly different scheme, according to whether T[i]<=bd. But what I wrote above should work.

Some other advice, if I may:

• Use `add` instead of `sum` in that routine, where you add up a finite number of values. It is faster.
• Instead of all those `subs` calls, you should be able to do it faster with `eval`. Now, I remember that you asked about it elsewhere, and that you do need the successive substitutions, in order. But why not substitute into the S[i-1] (to get "new" S[i-1] formulae) outside of the 'i' loop? That is to say, get your hands on fully explicit algebraic substitution formulae, newS[i],i=1..5, just once outside the loops. And then use just a single call to eval(X[i]*F||i(t),[..newS[i]=..] inside the loop. It would run much faster, I suspect. And the routine is slow.
• Remove N as a declared local and the first line which assigns to N using n, and change the parameter n::integer to N::integer:=1000.
• I would write procedures in 1D Maple input, not 2D Math. It looks bad in 2D. And I can actually see the delay which the GUI parses and displays the code. And it is wholly unclear visually whether the subscripted names are table references or not.
• Use the parameter names instead of args[i], all over, for consistency.

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

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