Why I am trying what I am trying in my ode program

This is the sort of thing that is very helpful. 

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 Fri, 8/10/12, Barton Willis <willisb at unk.edu> wrote:

> From: Barton Willis <willisb at unk.edu>
> Subject: RE: [Maxima] Why I am trying what I am trying in my ode program
> To: "Dennis Darland" <dennis.darland at yahoo.com>, "math maxima" <maxima at math.utexas.edu>
> Date: Friday, August 10, 2012, 6:29 AM
> > First there is probably a way to
> symbolically derive equations for the terms of Taylor series
> terms from the differential equations in terms of the
> initial conditions. 
> 
> For one first order DE, try something like:
> 
> (%i1) taylor_solve_de(de, x, y, [l]) := block([xo,yo,n,k],
> ? de : lhs(de)-rhs(de),
> ???xo : assoc('xo,l),
> ???yo : assoc('yo,l),
> ???n : assoc('n,l),
> ???taylor_solve_de_xx(de,y,x,yo,xo,1,n))$
> 
> (%i2) taylor_solve_de_xx(de,y,x,yo,xo,k,n) :=
> block([g,eq,sol],
> ???if k >= n then yo else (
> ? ? ???g : gensym("g"),
> ? ? ???eq : subst(y = yo + g *
> (x-xo)^k,de),
> ? ? ???eq : ev(eq,diff),?
> ? ???
> ? ? ???eq : tlimit(eq /
> (x-xo)^(k-1),x,xo),
> ? ? ???sol : algsys([eq],[g]),
> ? ? ???taylor_solve_de_xx(de, y, x,
> yo + subst(sol, (x-xo)^k*g),xo,k+1,n)))$
> 
> Taylor solution to y' = mu * y:
> 
> ? (%i3) taylor_solve_de('diff(y,x) - mu * y, x, y,
> xo=0,yo=1,n=4);
> ? (%o3) (mu^3*x^3)/6+(mu^2*x^2)/2+mu*x+1
> 
> Order nine solution to ' diff(y,x) = x*(1+y)^(1/2), x,?
> with initial condition (x=0,y=12)
> 
> ? (%i4) taylor_solve_de('diff(y,x) = x*(1+y)^(1/2), x,
> y, xo=0,yo=12,n=9);
> ? (%o4) x^4/16+(sqrt(13)*x^2)/2+12
> 
> What happened to those O(x^5) terms?? Let's try order
> 19:
> 
> ? (%i5) taylor_solve_de('diff(y,x) = x*(1+y)^(1/2), x,
> y, xo=0,yo=12,n=19);
> ? (%o5) x^4/16+(sqrt(13)*x^2)/2+12
> 
> Same approximation....oh, the DE has a solution that is a
> 4th degree polynomial :)
> 
> ? (%i6) subst(y = %, 'diff(y,x) - x*(1+y)^(1/2));
> ? (%o6)
> 'diff((x^4/16+(sqrt(13)*x^2)/2+12),x,1)-x*sqrt(x^4/16+(sqrt(13)*x^2)/2+13)
> 
> ? (%i7) %,diff;
> ? (%o7)
> -x*sqrt(x^4/16+(sqrt(13)*x^2)/2+13)+x^3/4+sqrt(13)*x
> 
> ? (%i8) scanmap('gfactor,%);
> ? (%o8) 0
> 
> (1) I thought I could generate the Horner form
> automatically, but I didn't get that to work.? I also
> think the code 
> ? ???could be modified to generate
> continued fraction solutions.? Variations are possible
> too (Frobenius-like solutions).
> 
> (2)? Of course, this could be boosted up to systems of
> first order DEs.
> 
> (3) This code will work poorly with floating point
> numbers---maybe the intent of your code was to do things
> like
> ? ???this, but stick purely to floating
> point calculations (for speed)?
> 
> (4) Likely my code is missing a few things that I can't
> think of right now.
> 
> --Barton
> 
> 
>