Timing comparison for diffeq program in c++, Maple and Maxima

In testing I've mostly stuck to problems which have analytic solutions. But most problems do not have analytic solutions at all. I think an analytic solution would always be preferable to a numeric one (I cannot think of any exceptions). But you cannot always get one. Sometimes, you cannot do the 'integrate' step. I've mostly been working on problems with analytic solutions for testing purposes. I do not question the superiority of maxima or maple for symbolic calculations. I have much more experience with maple. I also know there are other, more common, numeric methods. I got interested in this as a possibly overlooked one.

Dennis J. Darland
dennis.darland at yahoo.com
http://dennisdarland.com/
http://dennisdarland.com/dennisdaze/
http://dennisdarland.com/philosophy/
http://sode.sourceforge.net/
"According to the World Health Organization, the warming of the planet caused an additional 140,000 deaths in 2004, as compared with the number of deaths there would have been had average global temperatures remained as they were during the period 1961 to 1990. This means that climate change is already causing, every week, as many deaths as occurred in the terrorist attacks on September 11, 2001"
-- Peter Singer _Practical Ethics, Third Edition_, p. 216.


--- On Thu, 8/9/12, Rupert Swarbrick <rswarbrick at gmail.com> wrote:

> From: Rupert Swarbrick <rswarbrick at gmail.com>
> Subject: Re: [Maxima] Timing comparison for diffeq program in c++, Maple and Maxima
> To: maxima at math.utexas.edu
> Date: Thursday, August 9, 2012, 10:50 AM
> Dennis Darland <dennis.darland at yahoo.com>
> writes:
> > How, given the diffeq below, do you generate a Taylor
> series in 
> > Maxima?
> >
> > diff ( y , x , 1 ) = tan(sqrt(2.0*x + 3.0));
> 
> Well a 1d example like this is easy. Don't forget what a
> Taylor series
> is:
> 
> ? y(x) = y(c) + (x-c)*y'(c) + (x-c)^2/2*y''(c) + ...
> 
> You're saying that you have a formula for y'(x). In which
> case a
> (forward-Euler style) Taylor series is easy enough:
> 
> (%i1) rhs: tan(sqrt(2*x+3));
> (%o1)? ? ? ? ? ? ? ?
> ? ? ? ???tan(sqrt(2 x + 3))
> (%i2) assume(x>c);
> (%o2)? ? ? ? ? ? ? ?
> ? ? ? ? ? ?
> ???[x > c]
> (%i3) taylor ('integrate(rhs, x, c, x), x, c, 2);
> (%o3)/T/ tan(sqrt(2 c + 3)) (x - c)
> ? ? ? ? ? ? ? ?
> ? ? ? ? ? ?
> ???2? ? ? ? ? ?
> ? ? ? ? ? ? ? ?
> ? ? ? ? ? 2
> ? ? ? ? ? ? ? ?
> ? ? ? ? ???(sec (sqrt(2 c
> + 3)) sqrt(2 c + 3)) (x - c)
> ? ? ? ? ? ? ? ?
> ? ? ? ???+
> -------------------------------------------- + . . .
> ? ? ? ? ? ? ? ?
> ? ? ? ? ? ? ? ?
> ? ? ? ? ? ???4 c +
> 6
> 
> 
> (You need to quote integrate, otherwise Maxima tries to do
> the
> integration)
> 
> > or better, the system:
> > diff (x2,t,2) = 3.0 * diff(x2,t,1) - 2.0 * x2 -
> diff(x1,t,2) - diff (x1,t,1) + x1;
> > diff (x1,t,1) = 4.0 * x2 -? 2.0 * diff (x2,t ,1) -
> 2.0 * x1;
> 
> Well, if you look at the documentation for Taylor, you'll
> see that it
> can do multivariate Taylor series. You have a weird mixture
> of a second
> order and first order equation here. Presumably you want
> something that
> works for "any" equation you give it, so there's no point in
> trying to
> spot some structure in the example you've posted. As such,
> what you
> should probably do is transform this into a collection of
> coupled
> first-order ODEs.
> 
> Let
> 
> ? u1 = x1
> ? u2 = x2
> ? u3 = dx1/dt
> ? u4 = dx2/dt
> 
> Then you end up with
> 
> (%i1) eqs:
> ['diff(x2,t,2)=3*'diff(x2,t)-2*x2-'diff(x1,t,2)-'diff(x1,t)+x1,
> ? ? ? ? ? ? 'diff(x1,t)?
> =4*x2 - 2*'diff(x2,t)-2*x1];
> ? ? ? ? 2? ? ? ?
> ? ? ? ? ? ???2
> ? ? ???d x2?
> ???dx2? ? ? ? ? d
> x1???dx1? ?
> ???dx1? ? ???dx2
> (%o1) [---- = 3 --- - 2 x2 - ---- - --- + x1, --- = - 2 ---
> + 4 x2 - 2 x1]
> ? ? ? ???2? ? ?
> dt? ? ? ? ?
> ???2? ? dt? ? ?
> ? dt? ? ? ? dt
> ? ? ???dt? ? ?
> ? ? ? ? ? ? ? dt
> (%i2) subst(['diff(x2,t,2)='diff(u4,t),
> 'diff(x1,t,2)='diff(u3,t),
> ? ? ? ? ?
> ???'diff(x2,t)=u2, 'diff(x1,t)=u1],
> ? ? ? ? ? ???eqs);
> ? ? ? ???du4? ?
> ? ? ? ? ? ???du3
> (%o2)???[--- = - 2 x2 + x1 - --- + 3 u2 - u1,
> u1 = 4 x2 - 2 x1 - 2 u2]
> ? ? ? ???dt? ?
> ? ? ? ? ? ? ? dt
> 
> And the system is in fact linear here, so the Taylor series
> is trivial.
> 
> 
> I think that if you want Maxima to do something useful you
> should
> probably tell us what your "sode" code is actually supposed
> to do. Is it
> just a naive numerical integrator? In which case, give up on
> symbolic
> algebra packages and code it in C/C++ or Octave. You should
> also
> probably read some numerical analysis textbooks: just using
> Taylor
> expansions isn't the best approach.
> 
> Or maybe you want your code to do some cool maths that
> automatically
> spots special structures in the problems it's given? In
> which case, a
> higher level tool might be fantastic. But it really isn't
> clear what
> you're trying to achieve.
> 
> 
> Rupert
> 
> -----Inline Attachment Follows-----
> 
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>