Systems of Ordinary Differential Equations

by Dennis Darland Mail To: pal at dennisdarland dot com

Comments welcome

What omnisode is.
Where it came from.
Documentation.
Test Results.
Known Bugs.
Bugs Fixed.
Download.
Plans for the Future.
Items Accomplished.
Taylor Series - Bibliography.
Notes.
Old version of this page. (Before extensive revision on 12/21/2012)

Omnisode - What it is
1. It is mainly a Ruby program which generates a program to solve a set of one or more ordinary differential equations.
2. It currently can generate either Maple, Maxima, Ruby, c or c++ programs.
3. It uses a specifiable number of terms of the Taylor series of the equations.
4. The Taylor series is iteratively computed from recursive relations determined by the equations.
5. It runs under either Linux or cygwin (a linux like environment under windows).
6. Now it may be run under windows - but only generating Ruby.
7. A two shell scripts exist for each language for which a program can be generated.
8. Output for one script file goes to a file which is opened with emacs when the solution is completed. (This is all there is for Windows)
9. Output for one script file goes to a files which may be opened with a web browser after the solution is completed. (Not for Windows)
10. In the second case all the data may be viewed from a table.html.
11. A ".ode" file specifies any problem to be solved.
12. Note: where I say analytic solution, I really mean a solution given as a closed-form expression and evaluated at a point numerically.
13. It is open source software licensed under the GNU General Public License v3.0
14. Note: The peculiar way some of the code is written is due to the way arrays in maxima work (they seem to be call by value).
Omnisode - Where it came from
1. I became interested in this method from Y. F. Chang at the University of Nebraska, Lincoln in about 1977-1978.
2. Unfortunately, I became ill and was unable to complete my Masters degree.
3. I retained an interest in it, even though it was many years before I had a computer actually powerful enough for this method.
4. I intermittently remained in contact with Prof. Chang and Prof. George Corliss who had also been at the UNL.
5. Prof. Corliss is now at Marquette. I believe Prof. Chang is deceased.
Omnisode - Documentation
1. There is documentation on the ode files. (Revised 12/22/2012)
2. There is documentation on the script files used to run omnisode. (Revised 12/22/2012)
3. rdoc (Ruby Documentation) for release 0027 of the ode file preprocessor.
4. rdoc (Ruby Documentation) for release 0027 of omnisode. (diffeq generator)
5. rdoc (Ruby Documentation) for release 0027 of the diffeq post processor.
6. rdoc (Ruby Documentation) for release 0027 of diffeq.
Omnisode - Latest Test Results
1. Comprehensive Tests. (1/29/2013)
  The radius of convergence seems to be calculated accurately (which is all I use for now) but not the order as well especially for real singularities.
  I have not been able to figure out why, although I have spent quite a bit of time on it.
  Also the increment for a desired accuracy is not always calculated correctly. A maximum h may be set to help with this.
2. Tests of erf (Maple only)
  erf released (4/4/2013)
3. Tests of Si (Maple only)
  Si released (4/4/2013)
4. Comprehensive Tests. (4/24/2013)
  The radius of convergence seems to be calculated accurately (which is all I use for now) but not the order as well especially for real singularities.
  I figured out what is wrong, and don't know why radius works as well as it does. Order definitely gives poor results. See Notes. April 21, 2013
  Also the increment for a desired accuracy is not always calculated correctly. A maximum h may be set to help with this.
5. Tests using the ratio test to locate singularities (from Henrici Vol 1). (5/15/2013)
  This test was run on (and uploaded from) an old desktop that is very uncomfortable for me to sit at. I did not want to tie up the laptop I usually use.
  Thus I have not studied the results yet, though I did some testing in individual cases, Also excluded Maple and Maxima from the test, as these are not installed on the desktop.
6. Ppreparing for testing of location and order of singularities. Adjacent output in results file of info on radius and order for actual (if given in ode file), ratio test, three term test, and six term test. The results look promising. But order has trouble in three term test. Perhaps using more terms could help. I will make a release of this version.
Omnisode - Known Bugs
1. Sometimes computed increment (glob_h) does not give desired accuracy. (this may be impossible - I need to think more about it).
2. Comments in ode files start with "#" and must start line - instead of Maple "/* --------- */ syntax - which ode files are supposed to follow.
Omnisode - Bugs Fixed
1. First ode file in a batch is not copied to results directory. - Fixed.
2. Constant to full series power not accurate. (expt_c_sin_new.ode) - Fixed.
3. Linear function to full series power not accurate. (expt_lin_sin_new.ode) OK except h determined is too large for desired accuracy.
4. Full series function to full series power not accurate. (expt_sin_sin_new.ode) OK except h determined is too large for desired accuracy.
5. Fixed missing "tail.html" from release. (in rev0027)
Omnisode - Download
1. The current omnisode (0035) was published on the web May 20, 2013.
2. Both tar.gz and zip versions - contents should be identical.
3. The latest published version may be obtained at The Sourceforge Download page for Omnisode
4. Omnisode download stats as of 4/15/2013
Omnisode - Plans for the future
1. It seems in some cases spaces are needed to separate tokens in equations in ode file.
2. Need better error messages on syntax errors is equations.
3. Individual equations in ode file (1st section) must be on one line. Should permit multiple lines.
4. Investigate use of interval arithmetic. Now (12/15/2012) think this is probably impractical.
5. Ruby support using APFP Almost certainly impractical.
6. Perl support - Little reason to do - except I like playing with languages.
7. Python support - Little reason to do - except I like playing with languages.
8. Note that the support would be for the programs generated, omnisode will remain a Ruby program.
Omnisode - Items accomplished
1. May 22, 2012 - Language independent ".ode" files (they define a set of equations to solve)
2. May 22, 2012 - Documentation of ode file format
3. May 24, 2012 - Capability to add test data to a html table automatically
4. May 25-31, 2012 - Links to detail files added to tables generated.
5. May 25-31, 2012 - Automatic generation of test scripts.
6. May 25-31, 2012 - Automatic indentation of programs generated.
7. May 25-31, 2012 - Work files moved to "temp" subdirectory.
8. May 25-31, 2012 - Added "ode.over" file to override values in all ".ode" files for batch tests.
9. June 5-12, 2012 - Improved significantly handling of systems of (multiple) equations.
10. June 12, 2012 - More documentation.
11. June 15-18, 2012 - Optimized factorial calculations.
12. July 2-10, 2012 Treat constants and linear expressions as special cases of series in order to optimize these common cases. More efficient bothe in time and space. (Mostly Done - except difficulties in a few cases).
13. July 11, 2012 Put in warnings for special cases which seem to have problems - or have trouble checking. (Also need to check these cases again).
14. August 13, 2012 C++ code generation functioning (Almost c also).
15. August 21, 2012 Now C as well as C++.
16. September 6, 2012 Made to support negative increment (H).
17. September 6, 2012 Added capability to display values calculated only at specified intervals.
18. November 7, 2012 Now Ruby code can be generated.
19. December 14, 2012 Attempting to control error using estimate of radius of convergence and other factors - still needs much testing.
20. December 22, 2012 Revised documentation.
21. January 16, 2013 Worked on special cases with difficulties - devised tests for which I had analytic solutions and fixed any problems.
22. January 26, 2013 Examined how estimation of singularities is working. Radius works well. Order less accurate.
23. January 27, 2013 Added link to ode.over above rest of data in output table.
24. January 27, 2013 Devised way to set maximum h smaller than display_interval.
25. January 27, 2013 Cleaned up code (Mainly deleting dead (commented out) code).
26. January 27, 2013 Cleaned up ode files (Mainly removing obsolete glob_h)
27. February 24, 2013 Added support for Windows - just for Ruby.
28. April 4, 2013 Support for Si and erf functions (Maple only) released.
29. May 15, 2013 Added logic for ratio test for singularities.
30. May 20, 2013 Results file now contains adjacent output data of singularities - actual (if location given in ode file), ratio test, three term test, and six term test.
Taylor Series - Bibliography
1. D. Barton, I. M. Willers, and R. V. M. - "The automatic solution of ordinary differential equations by the method of Taylor series", in The Computer Journal, 14(3): 243-248, (1970)
2. Peter Henrici - Applied and Computational Complex Analysis: Volume 1: Power Series-Integration-Conformal Mapping-Location of Zeros(1973)
3. Y. F. Chang - "Automatic Solution of Differential Equations" in Constructive and Computational Methods for Differential and Integral Equations, Springer-Verlag, (1974)
4. Peter Henrici - Applied and Computational Complex Analysis: Volume 2: Special Functions-Integral Transforms-Asymptotics-Continued Fractions(1977)
5. George F. Corliss - "Integrating ODE's Along Paths on Riemann Surfaces", in Proceedings of the Seventh Manitoba Conference on Numerical Mathematics and Computing (September 29 - October 1, 1977)
6. George Corliss and David Lowery - "Choosing a stepsize for Taylor series methods for solving ODE's" in Journal of Computational and Applied Mathematics, volume 3 and 4, (1977)
7. Y. F. Chang - Start of Draft of book on Taylor Series (1978)
8. D. Barton - "On Taylor Series and stiff equations" - ACM Trans. Math. Softw. 6, 3 (Sept. 1980).
9. Louis B. Rall - Automatic Differentiation: Techniques and Applications, Springer-Verlag, (1981)
  (this was hard to get, I finally found a discarded copy from the author's university's library on a web bookstore)
10. George F, Corliss - "Integrating ODE's in the Complex Plane - Pole Vaulting", in Mathematics of Computation, Volume 15, Number 152, (October 1980)
11. Y. F. Chang and G. Corliss - "Ratio-Like and Recurrence Relation Tests for Convergence of Series", in J. Inst. Maths Applics 25, (1980)
12. George Corliss and Y. F. Chang - "Solving Ordinary Differential Equations Using Taylor Series" in ACM Transactions on Mathematical Software Volume 8 Number 2, June 1982, pp. 114-144
13. Y. F. Chang - "The ATOMCC Toolbox" in BYTE, Volume 11, Number 4, pp 215-224, April, 1986
14. Peter Henrici - Applied and Computational Complex Analysis: Volume 3: Discrete Fourier Analysis-Cauchy Integrals-Construction of Conformal Maps-Univalent Functions(1986)
15. Andreas Griewank and George F. Corliss - Editors - Automatic Differentiation of Algorithms: Theory, Implementation and Application, siam (1992)
16. Y. F. Chang - ATOMFT User Manual, ver. 3.11 (1994)
17. Martin Berz, Christian Bischof, George Corliss and Andreas Griewank - Editors - Computational Differentiation: Techniques, Applications and Tools, siam (1996)
18. Angel Jorba and Maorong Zou Taylor User's Manual, November 13, 2001
19. George Corliss, Christele Faure, Andreas Griewank, Laurent Hascoet and Uwe Naumann - Editors, Automatic Differentiation of Algorithms: From Simulation to Optimization , Springer, (2002)
20. Angel Jorba and Maorong Zou A software package for the numerical integration of ODEs by means of high-order Taylor methods, January 27, 2004
21. Martin Bucker, George Corliss, Paul Hovland, Uwe Naumann and Boyana Norris - Editors Automatic Differentiation: Applications, Theory, and Implementations, Springer, (2006)

Notes

May 21, 2013
I ran general tests yesterday and today. It was a comedy of errors. I had changed the path to the results, from html/results to html/omniresults, in order to better specify that these file need not be backed up. I thought I had changed it everywhere. But I hadn't in the 1st line of gentest.sh, which removes the prior test. Then all the results were not saved. Then I corrected that and reran it. This time, in the paths in the t*l.sh files I had accidentally changed it in the results file name itself. I wrote a quick ruby program and renamed all the files. Then I looked at the results. C and C++ were giving inaccurate results on add_c_lin.ode, but Ruby, Maple and Maxima gave good results. I am too tired to work more on it now. I remember I haven't changed the path in the ruby script for windows either.
Later - Have new laptop set up ruby works OK - but duplicated problem with c and c++ - maxima and maple not installed yet.
May 20, 2013
More work on preparing for testing of location and order of singularities. The results look promising. But order has trouble in three term test. Perhaps using more terms could help. I made a release of this version at Download omnisode (0035).
May 19-20, 2013
Modified omnisode to accept (fully implemented for time being only for generating ruby - will need to preprocess subscripts) location of closest singularity and compare actual (given) radius and order with those computed by the ratio, three term and six term tests. These are output in the results files to which there are links in the results table
Now that these results are easily examined, I can start some serious debugging!!!
May 19, 2013
I am getting ideas on what to do next. I am considering adding columns to the results table for the actual radius and order of singularities, and also the results of the equations derived by Corliss and Chang. The radii given by those seemed close, even though I think I found an error in their derivation. (also retaining the radii given by the ratio test.)
May 15-16, 2013
Worked on implementing ratio test (see Henrici Vol 1) instead of Corliss' and Chang's tests. See Test Results..
May 14-15, 2013
Worked on debugging ratio test.
Main difficulty was I was only keeping minimum radius - instead of radius for each series.
I was making progress, but became exhausted and went to bed.
May 13-14, 2013
Modified omnisode to use ratio test as found in Henrici, Volume 1
Debugged and started comprehensive test (on older slower computer - I need this one for other things)
NOTE: Also found I was getting bus error with emacs - switched to emacs-x11 and resolved the issue.
Very Tired.
May 3, 2013
Version 0033, matching latest test results, was published 4/30/2013.
This page was not updated till today due to technical difficulties.
April 29, 2013
Someone has recommended Chebyshev polynomials.
Some links are:
1. Chebfun
2. One hour tutoriial(I couldn't get to work)
3. History
4. chebfun
5. Approximation Theory and Approximation Practice
I ran the problem given in the 1st link with omnisode - Test from Chebfun
It has trouble with expt near zero, but does OK, it seems, with the products.
The link gave an exact symbolic solution but not the numerical value.
April 24, 2013
Looked over test results and uploaded them. See Test Results.
Cases intended for testing singularities (singN.ode (where N is a digit) and nonlinearN.ode (same) gave accurate results for radius.
I did not get good results for order and it now says "Unknown" - that was easier than removing it altogether.
Some functions which I would have expected to be entire gave finite positive radii.
The equation will always give an answer and my only check is that the radius is positive.
I had been checking for positive order, but all orders except negative integers are permissible, and how can you check that with floating point?

After more study:
I kept getting wrong answers for the order, though the radii were close in my test cases. I kept studying the analysis in Chang's typescript more & more closely.
The trouble is going from equation [11] to equation [12]

d f(x)/dx = - p * (u(x)**(-p-1)) * d u(x)/dx [11]
to
F(2)*U(1) / h = - p * F(1) * U(2) / h [12]

I think there is an error in that one needs to evaluate the derivatives of u(x)**(-p-1) at x=0. and I think Corliss took the derivatives of u(x) at x= 0 and took them to the (-p-1) power. Though I cannot quite get his result. (I have not tried very hard as it would have to be wrong). But if one takes the derivatives of [11], then the terms of u(x)**(-p-1) will not disappear as they do in [13], unless p is a negative integer.
[Added 4/25/2013 - the equations (on which the "three term analysis" (p. 79.)) referred to are in part 4 (pp. 74-75) of Chang's typescript.
Also it should be noted that the "six term test" (p. 96) is based on the relations used in the three term test (p. 82).
A variation of the results also were published in Y. F. Chang and G. Corliss - "Ratio-Like and Recurrence Relation Tests for Convergence of Series", in J. Inst. Maths Applics 25, (1980), though not all of the derivation.]
April 23, 2013
Checked over test results.
Reduced glob_max_h for a few cases.
Removed mtest6 cases with larger glob_max_h
Adjusted bounds for a few cases.
Started another test.
April 22, 2013
Studied results of yesterday's test.
Hard to check results, because with different speeds of different languages, they timeout at different points.
Modified program to report data at point where singularity is first detected in order to have data in table easier to compare.
Started another test.
April 21, 2013
Went over the derivation of equations for poles.
I believe there is an error going from eq [10] to eq [12] in Chang's typescript (see bibliography - it's in part 4 ).
There eq [12] is:
F(2) * U(1) / h = - p * F(1) * U(2) /h
I believe it should be:
F(2) * U(1) / h = - p * F(1) * (U(2) ** (- p - 1)) / h
The derived equation for the radius of convergence seems to work anyway.
I haven't figured how to derive results from the corrected equation [12].
I wasn't using the order anyway.
I started a new comprehensive test. The order is now given as "Unknown".
April 19, 2013
Worked much trying to get radius to work for real poles.
Redid formulae and couldn't make it work any way I tried.
Worked many hours while following news on search for Boston Marathon bomber.
Did however make some fixes, though they did not fix basic problem.
Not enough for a release.
April 18, 2013
Worked on calculations of poles.
Made some corrections.
Complex poles - both radius and order working very well.
Real poles - radius works very well - but not order.
Suspect problem may not be in calculation itself - need to check it - and use of arrays - etc.
April 4, 2013
rev0030 released - Si and erf (Maple only)
March 30, 2013
Added and tested Si (Maple only)
March 29, 2013
Finished and tested erf (Maple only)
March 28, 2013
Added logic for erf for constant and linear arguments. (Will only work for Maple - needs erf built-in function.)
Started testing it.
Still need full series case - shell added.
Also have shell for Si function added. (Also for Maple only - requires Si built-in function.)
March 19, 2013
Studied some generated code and did some testing. (All OK)
Thinking about how to go about handling singularities.
March 2, 2013
Decided to focus next on testing and checking logic for locating singularities.
February 25, 2013
Created zip archive and added to sode sourceforge download section.
Contents should be identical with tar.gz file.
Added for convenience of windows users.
February 24, 2013
Finish changes to also support windows (rev0029) - but just Ruby code generation.
I don't have the other languages under Windows.
February 23, 2013
Modified omnisode.rb to use "File.join"
Had difficulty with generated code so "entry.html" is now generated in current directory."
Scripts will have to be modified for this.
Later:
Have ruby program tr.rb working under windows (functions as tr.sh under Linux and cygwin)
Will have to modify (and test) shell scripts to account for entry.html being moved.
Assumes editor in path named "e" - my name for emacs under windows.
February 22, 2013
Started working on making omnisode work under windows.
I'm using windows 7 and installed the latest Ruby today.
Installed omnisode0027 on windows machine.
First just working on generating ruby. (I don't have any other languages installed on windows.)
Only change to code so far is path separator (which would also work with other languages).
Used 'SEPARATOR' to indicate '\' where linux (& cygwin) would use '/'.
Haven't decided how program should determine how to set that variable.
Made bare minimum batch file to test with.
Omnisode worked with sin.ode and ode.over.
Omnisode failed to solve Sin.ode did not work right without ode.over. (Needs investigation - probable happens under linux also.)
Need substitute for linux shell scripts (Ruby program?)
February 20, 2013
Added rdoc (ruby documentation) for the various programs in used in omnisode.
February 3, 2013
Test downloaded latest version of omnisode on a different machine (which uses cygwin).
Only difficulty was emacs is not set up same on this machine.
Only tested with Ruby, but "ruby diffeq.rb" worked after trying "tr.sh" which had created diffeq.rb.
I don't have any of the other languages on that machine.
January 30, 2013
The code seems stable, and I need to attend to other things.
I will let it incubate in the back of my mind before I do more - this is how I get my best ideas.
January 28-29, 2013
Did testing and uploaded table of results.
Comprehensive Tests.
The radius of convergence seem to be calculated accurately (which is all I use for now) but not the order as well especially for real singularities.
I have not been able to figure out why, although I have spent quite a bit of time on it.
Also the increment for a desired accuracy is not always calculated correctly. A maximum h may be set to help with this.
Made release rev0027 of omnisode 1/29/2013.
January 27, 2013
Added link to copy of ode.over to output table.
Made first ode file copy correctly.
Modified singularity code so it tells which equation.
Made so glob_max_h can set h less than would be otherwise computed - and which may be less than display_interval.
Cleaned up both omnisode.rb and ode files.
Started comprehensive test in all 5 languages.
January 26, 2013
Ended up using equation from Chang and Corliss (1980)
Results added to test results section.
January 25, 2013
Attempted to fix problem with order of pole for real singularities.
No luck, although I did get re-familiarized with the code.
Spent about three hours.
Later: Worked on this.
Studied - Y. F. Chang and G. Corliss - "Ratio-Like and Recurrence Relation Tests for Convergence of Series", in J. Inst. Maths Applics 25, (1980)
Derived related equation (using Maple). Different because they used different convention for indices. Better results but still not very good.
January 24, 2013
Studied test results for radius of convergence and order of poles.
Both appeared correct for complex poles.
Radius of convergence but not order of pole seemed correct for real poles.
I could not locate an error in the program.
January 17, 2013
Released rev0025 of omnisode.
January 15-16, 2013
Devised tests for cases I didn't have analytic solutions for comparison. Five of the eight worked well. The other three need attention.
New Tests for eight special cases.
Later: Initially 3 had problems - due to problems in entering the equations in the ode files. Two required a smaller h for better accuracy.
New Tests for seven special cases.
Two cases requiring smaller increment for good accuracy.
January 14, 2013
A Summary Table (1/14/2013) With a few bugs fixed - Extensive tests for single equations using c. c++, Ruby, Maple, and Maxima.
Single equations (not yet systems) work fairly well.
Tests with multiple equations (1/14/2013) still one problem is less accurate than others.
These tests just using Ruby - to save time.
Tests for multiple equations using Ruby.
The current omnisode (0024) was published on the web January 14, 2013.
The latest published version may be obtained at The Sourceforge Download page for Omnisode
January 12, 2013
Had time to take a quick look and found the problem in c.
It was introduced by one of the last of the other fixes.
I corrected it and started another big test.
Later: Looked over test results carefully.
All lines in test script resulted in rows in results table.
mult_sin_lin.ode was not as accurate as others.
Found mistake and started another test.
January 11, 2013
Test finished - some sort of problem with script (it seems) for c.
Too busy today to work on it.
January 10, 2013
Fixed remaining 3 problems. Actually just 2 hard to find bugs.
Put back all commented out tests in script and started test on single equations.
Need to look closely at results - especially results on singularities.
January 9, 2013
I had previously just checked for test results in the output tables mainly for low numbers of correct digits.
I had noticed some results missing - I knew these could be caused my syntax errors in the generated code - but did not think there were many.
Also there were a few (just for Maxima) which I had removed from the script which would get "stuck."
Today I went though and noted all rows for which there were some results, but some rows were missing.
I found 49 rows missing.
I am going to start making corrections tonight, but it will probably take a few days to make corrections, run tests, and make a new release.
After working through the evening, I had fixed all except 3 cases of missing rows. Required only 4 fixes I think.
There is also, though, a problem with estimation of singularities. I may have introduced it in one of the 4 fixes.
That area has been less tested anyway. Everything else has to work right before it can.
January 6, 2013
Uploaded release 0023
January 5, 2013
Uploaded results with fix from Jan 3 uploaded.
Will make release with fix soon.
January 3, 2013
Found and fixed problem in emit_pre_assign and emit_assign.
Re-running tests.
December 22, 2012
Uploaded rev0022 of omnisode.
Also extensive test results on single equations uploaded.
December 21, 2012
Test equation is:
y''''' = -y'''
Error appears in y''''
Somehow must miss its calculation or using some uninitialized value in its calculation.
Will sleep on it!
NO SLEEP!!! - Worked ALL Night!!!
I think I fixed the problems I had detected!
8:30am Started General Test.
Will sleep (hopefully) while test runs.
Single equations seem to work fairly well.
Some problems with multiple equations.
December 20, 2012
I studied and improved debug code.
Located and fixed errors in pre_emit_assign and emit_assign.
Can tell these are correct 1st time through.
There must be additional error in jump for higher order equations.
Began inserting debugging there.
December 19, 2012
I decided to do testing on a problem that had higher order derivatives on both lhs and rhs, but a single equation.
Decided on:
y''''' = - y'''
I put in debug to give both calculated and evaluated values of analytic solutions of first several terms.
At the beginning of a run, they were very accurate, but then suddenly a term goes wildly astray.
Thought it might be hardware.
Tried on a different computer - had to install ruby and program.
Had similar results - but not exactly same. - Sudden deviation - but not same place.
Must be software bug - like memory clobber - but some of the languages I am using are supposed to check that.
Enough for now!!!
December 17, 2012
Tests with different display interval for a problem with multiple equations.
However computed glob_h for desired accuracy was generally smaller than interval anyway.
See Special tests 0002
December 14-15, 2012
Completed general test and one special test comparing different interval sizes for the same problem.
See General test 0001
In the following table results are given for various increment sizes.
The increment is no longer set directly - the value used is computed to try to achieve a desired accuracy.
But the increment will not be larger than an interval set to display the results at. So there is still some control.
It turns out larger increments are better - at least up to some point. I still need a lot more testing on this.
See Special tests 0001
Note: Results files larger than 10MB have been removed. Many are still larger than 1MB. Some had been over 500MB.
December 13, 2012 (Later)
1. Fixed time in c and c++.
2. Discovered larger error when using larger interval seems more due to larger range tested than increased step error.
3. Started work on testing specific cases for c only (ones showing larger error in general test).
4. Rerunning general test with time fixed in c and c++.
December 12-13, 2012
Looked at test results and ran a few more smaller tests.
Observations:
1. Because some tests reached max time, Maxima has greater observed accuracy, but only because it is slower and had fewer iterations to accumulate error.
2. Times in c and c++ are too large. I can tell the tests are running much faster than measured time. Maybe what I am using as a second is really a millisecond?
3. Most tests run accurately.
4. Most tests that are not accurate, are accurate when smaller increment (h) is used.
5. When the smaller h is used, the results seem too accurate to be mistaken.
6. Need more systematic testing for cases where error seems larger (with similar h).
7. Note: because results often time out, before completing range, increasing h not only increases error on individual step, but also greatly increased the total range covered - thus magnifying the error reported.
8. Regarding this last note: the last fact may be overcome with more systematic settings in ".ode" and "ode.over" files.
December 11, 2012
Looked at problem again with correct debugging.
It appeared that I had already solved problem,
and that remaining appearance of problem was in error in debug code.
Started another general test for overnight.
December 10, 2012
Fixed code added for debugging.
Then too tired to debug problem.
Also updated project for SourceForge.
December 9, 2012
Isolated one error, but could not figure it out - could not seem able to think clearly.
December 4, 2012
Fixed a couple things related to problem below.
Seems I have, at least, a problem with a sign somewhere.
November 28, 2012
Located one problem and fixed.
Emit_pre_diff needed modification for the special case handling.
Still something wrong.
Too tired to think any more.
November 27, 2012
Narrowed bug down.
I believe it is in emit_pre_assign and emit_assign (at least).
But I could not think clearly about how it should work.
November 22, 2012
Added debug code for 1 specific problem (diff2.ode)
Debug code seems OK, but still haven't solved problem.
Seem unable to think (clearly)
November 21, 2012
Added Total Expected time to solve problem to output and replaced expected lime left with it in table produced.
I think it is more useful.
November 19, 2012
Seems there definitely is a problem (see 11/18/2012)
Studied code and didn't find the problem.
May have to resurrect some debugging code.
Too tired tonight.
November 18, 2012
I went over test results.
It seems that there are accuracy problems iff the derivative on the lhs is order greater than 1.
That is an isolated section of code to examine, but I am too tired today.
November 16, 2012 (early morning)
Fixed problem in ode files with lines containing only blanks - introduced with logic for comments.
Improved logic for determining size of increment - but still not really correct.
Started another general test.
(later) Went over results, fixed some problems and started retest.
November 11, 2012
Comments now OK in ode files (except for top equations section)
Increment automatically calculated to obtain desired accuracy for specified range. (Including negative increments.)
The last item is coded and runs but it seems not to always give the desired accuracy and thus needs much more work and testing.
November 10, 2012
Continued working on the same things - some progress
November 9, 2012
Worked mostly on automatically determining optimal increment size.
It started working, but then I got sidetracked on getting comments working correctly in preprocessor (preodein.rb).
Making progress there also.
November 8, 2012
Working on testing.
Tweaking ode files and not using "ode.over" - thorough testing requires that for now. (Especially size of increment.)
Would be good to make program select values itself - hopefully that will come later.
November 7, 2012
Worked on logic to adjust step size based on estimate of radius of convergence.
Started test on this including c, c++, Maple, Maxima and Ruby.
November 6, 2012
Indenting now OK.
Started lengthy general test.
Intend to make some support to adjust increment size for singularities.
Also to take larger step size as long as sufficiently accurate.
Also to support Ruby with my APFP library.
November 5, 2012
Have Ruby code generation working.
Also working on making output code indented properly for every supported language. (Maple, maxima, c , c++, Ruby).
November 4, 2012
Studying articles relating to increment size and pole-vaulting around singularities.
Interesting - but seems (especially pole-vaulting) to require special analysis of individual cases.
Started adding support for generating Ruby code.
October 23, 2012
The tests with the larger increment (below) have catastrophic errors if they pass over singularities. (I have known this,
but failed to comment on it here.) I have yet to work out how to handle this - I am researching and thinking about it.
September 21, 2012 Tests with larger increment (0.005) - Most results OK - But some very bad - No changes in code.
Report of tests for c. c++, Maple, and Maxima (WITH LARGER H - 0.005).
Expect to mostly be studying problem of error control and experimenting with code - no more releases likely for a while.
September 8, 2012
The current omnisode (0018) was published on the web.
September 7, 2012
Improved September 6th work and made test below:
Tests with negative H and/or display interval.
September 6, 2012
Made to support negative increment (H)
Added capability to display values calculated only at specified intervals.
September 3, 2012
Improved correct digits reporting (had been incorrectly reporting value for last equation for all in systems of equations)
Better report of tests for c. c++, Maple, and Maxima.
Released this version (0017).
This version (0017) may be obtained at The Sourceforge Download page for Omnisode
September 2, 2012
Fixed problems in my work done on August 30, 2012. Started general test.
August 30. 2012
Fixed problem with "Correct Digits" (see below) and started general test.
August 22, 2012
I have realized the "Correct digits" reported in the cases of multiple equations is that of the last equation - for all of the equations.
This will be fixed in the next release.
August 21, 2012
Ran general batch test with fixes from yesterday
There were some problems, so there is more to do before I post results.
However I can give better time comparisons.
I discovered the time reported by Maple is CPU time - so I am using the difference in the times in the time stamps instead.
Maxima can report either clock or CPU time - there are 2 functions - I am using clock time.
C and C++ report clock time.
The following was obtained with mtest6.ode
Maple did 100 iterations in 24 seconds (17 seconds CPU)
C++ did 69 iterations in 60 seconds
Maxima did 8 iterations in 67 seconds
This is:
Maple 4.17 iterations/second
C++ 1.15 iterations/second
Maxima 0.12 iterations/second
NOTE: The comparison is not altogether fair. Digits was set to 32 for Maple, while only 16 digit precision was obtained for C++ and Maxima
After Sleep:
Fixed problems with changes to support local variables in ode files for maxima.
Fixed Digits reported in c and c++ (got it right this time).
Reversed order of entries is results table.
Added script files for c.
A Summary Table (8/21/2012) Now including c and c++ (really just c code but compiled with c++) is at:
Extensive tests for c. c++, Maple, and Maxima.
Published current omnisode (0016) on the web.
The latest published version may be obtained at The Sourceforge Download page for Omnisode
August 20, 2012
Fixed a number of problems with last release (omnisode0014)
1. Made to work with c as well as c++ (using gnu gcc and g++)
2. Fixed time measurement in c and c++ (Also note Maple time is CPU time - Maxima, C and C++ time are clock time)
3. Fixed digits field in table for c and c++
4. Added number of correct digits to output and table when analytic solution given
5. Added beautify program for tcppl.sh script.
6. Made preodein.rb handle local declarations for maxima (only affects analytic solutions)
August 12-13, 2012
Testing and debugging c++. Test results at Extensive tests for Maple, Maxima, and c++
For some reason the tests "expt_c_sin", "expt_lin_sin", "expt_sin_lin", and "expt_sin_sin" got stuck in Maxima ands were removed from the test.
Need to study results thoroughly.
August 11, 2012
Ordered the 3 volumes of Peter Henrici - Applied and Computational Complex Analysis
Finished work on generating c and c++
Started testing c and c++ - a few corrections, but going well.
August 10, 2012
Added blocks and declared all local variables in maximas. Maxima was a little slower with the blocks added, but not very significantly. I knew where they all went, as I had them (local declarations) for Maple, and mint checks that so I should not have missed any. I have found out how to use big floats in maxima. It will be quite a bit of work, and as maxima is sort of slow, and presumably this would make it slower, I am not pursuing this now. Most of the time, I suspect input data to real problems is not known that accurately anyway. I will be looking into interval arithmetic though - probably not for maxima. Maple has something but I don't know if it will do what I need. I understand what mode_declare does now, but it doesn't support big floats, so I will avoid it for now. I don't have not had problems without it, and I think it is only effective if the code is compiled. I am not interested in compiling as then I could not use big floats, and I want to do that someday. Next, I guess, is finishing up C and C++.
August 9, 20121
I ran a timing test comparison for c++, Maple and Maxima.
Set time limit to 1 minute.
Maximum iteration was sufficient to solve entire problem.
Problem was sin.ode for all three.
Output was redirected to file only - none to screen
Results:
c++ solved entire problem in 23 seconds (140001 iterations)
Maple performed 11542 iterations in 60 seconds
Maxima performed 372 iterations in 60 seconds
That amounts to:
c++ 6,087 iterations/second
Maple 192 iterations/second
Maxima 6.2 iterations/second
c++ is almost 32 times faster than Maple which is almost 31 times faster than Maxima.
So c++ is almost 982 times faster than Maxima.
Later (I think my dates may be off the last few days):
Removed mode_declare's from maxima
Then maxima performed 394 iterations in 60 seconds
Or:
6.6 iterations/second.
Was able to get bloats (with trig functions) to work in maxima. I had misunderstood syntax before.
Am going to try to compile maxima code - I realize I cannot do this with if I use bfloats - which I don't now.
August 8, 2012
c++ mostly working. So much faster will have to modify output.
Files produced are too large to view, if enough iterations are used in order measure time.
August 7, 2012
c++ version running and giving good answers with very limited testing.
Some problem with displaying time data, but it appears to be much faster than either Maple or Maxima.
August 3-6, 2012
Working on adding c and c++ to languages generated.
Going much easier than anticipated.
Mostly just tedious parts left - a few tricky ones also.
July 27, 2012
Checked coding of Chang's case 5 for expt_lin_lin.ode and found nothing wrong. Tried using Taylor series of 250 terms but this did not help. The analytic solution and derivation of Chang statements are beyond my skill level currently. I got to thinking that I remembered something on this in Rall (1981). Looked and found it. Not sure if it will help or not. No time now.
July 11, 2012
Made release 0013 It may be downloaded from The Sourceforge Download page for Omnisode- Main change is to handling constants and linear expressions as special cases.
See test results above
July 1, 2012
Did sqrt, and modified array allocation to save memory as made possible with optimized logic.
June 27-30, 2012
All cases of exponentiation coded, except sqrt.
Next - a lot of testing!
June 25, 2012
Continued changes for constants and linear expressions - exp, and log done.
Log has truncation problem

June 24, 2012

Language	Compute every time	Pre-compute all	Compute as needed and save	With Linear Opt
Maple	1 min 40 sec (16 digits)	29 sec (16 digits)	41 sec (32 digits)	20 sec (32 digits)
Maxima	40 min 33 sec (16 digits)	31 min 58 sec (16 digits)	31 min 25 sec (16 digits)	8 min 52 sec (16 digits)

What I used here was the time left - I should have used the time elapsed + time left. - I did not think of that before.
But I don't think it really matters - they should be proportional to each other.
By linear opt, I mean taking advantage of the fact that the operand to the sin function is linear (In this particular example).
Did more testing on arccos & arctan. Fixed 2 typos in arctan. Both did better when dividing increment by 100 even though multiplying iterations by 100. Thus it must be truncation error. (Getting about 20 place accuracy.)

June 23, 2012
I ran into less accuracy with arccos and arctan than I have been typically getting. (I know the domain of values matters and that I need to do much more thorough testing. I am to mostly trying to catch blatant errors now. But I was getting significantly less accuracy (for the domain I was using) for arccos and arctan than other functions. So I decided to look at the Fortran code generated by Chang's Atomft. First I tried it on Linux, but was unable to find g77. So I tried cygwin, and after a little work (mainly running d2u on the files and reworking the make files to my liking. - I had used Atomft before but only retained the original zip file.) ran it on the arc functions. I got error messages - they are not supported by Atomft. Then I went back to Linux and, with a little web browsing, found out that g77 is now called gfortran - and got the same results.
June 21, 2012
Continued changes for constants and linear expressions - sinh, cosh, tanh, and arcsin done.
June 20, 2012
Continued changes for constants and linear expressions - division, cosine and tangent done.
June 19, 2012
Continued changes for constants and linear expressions - subtraction and multiplication done.
Also looking at some Taylor series capabilities in Maxima - should also look at Maple
June 18, 2012
Same tests, but factorials computed only as needed, but then saved for when needed again should be same results except faster(June 18, 2012) Started revisions to treat constants and linear expressions as special cases of series in order to optimize these common cases.
I have modified the basic structure for this and implemented and tested for addition and sine.

June 17, 2012
Same tests, but factorials pre-computed should be same results except faster(June 17, 2012)
I have revised the factorial logic again. Instead of pre-computing all the factorials, I am just computing them the first time needed. Also I noticed (Affecting maple only) I had Digits set to 16 in sin.ode, to 50 in h2sin.ode, h3sin.ode and sing4.ode, and to 100 in sing4.ode.
These are all now set to 32. (like everything else).
The timing results for sin.ode were:

Language	Compute every time	Pre-compute all	Compute as needed and save
Maple	1 min 40 sec (16 digits)	29 sec (16 digits)	41 sec (32 digits)
Maxima	40 min 33 sec (16 digits)	31 min 58 sec (16 digits)	31 min 25 sec (16 digits)

I have another general test running. It will take several hours and is on another computer so my editing here does not affect it

June 16, 2012
Same tests, but just increment (H) increased another factor of 10 new results (increment increased by factor of 10 )(June 16, 2012)
NOTE: sometimes maxima's error is less - but because it reached the maximum time (set to 15 minutes) before it reached maximum iterations.
Pre-calculated factorials to increase efficiency. Running test.
June 13, 2012
Same tests, but increment and number of iterations both increased by factors of 10 new results (mult H & and iterations by factors of 10 )(June 13, 2012)
Same tests, but just number of iterations increased another factor of 10 new results (iterations increased by factor of 10 )(June 13, 2012)
June 12, 2012
After much thought and work on systems of equations I started getting better results.
The poorest are mtest5.ode, mtest9.ode, mtest9_rev.ode and log.ode.
BTW: mtest6.ode and mtest6_rev.ode were previously called complicated.ode and complicatedrev.ode
June 8, 2012
I added a capability to give initial conditions of order for a dependent variable higher than what appears on the LHS, to go along with permitting them on the RHS. This, I think is needed for some systems of equations. For such terms, there values will not be updated in the "atomall" loop. Rather, they will be projected in the "jump_series" function. I have all this coded and running, but have not yet made any examples where it is really used. I had to make a change to the "ode" files in order to accomplish this. Next I need to update the existing ode files (I have only done one.), and come up with some real examples.
June 6, 2012
Made global variable in omnisode to set which of 2 ways to jump - relating to the subiter methods and multiple equations. It runs but too tired to really test. Need better sets of equations - I have some ideas - but too tired now.
June 5, 2012
Had to make changes to omnisode to handle the sort of problem I was thinking of - I think I need a yet harder problem!
Table comparing tests using three methods on difficult system of equations
I realized I made a mistake on the prior table - I had the method overridden in "ode.over." I reran it and got these results:
Table comparing tests (2nd try) using three methods on difficult system of equations
I need interdependent equations at higher orders I think. I think maybe I have now made it so that one need not have the highest derivatives on the LHS!
June 4, 2012
Uploaded omnisode0010.tar.gz to The Sourceforge Download page for Omnisode Modified omnisode to use either one sub-iteration or total order sub-iterations or total order + max terms sub-iterations according to the variable glob_subiter_method is set to 1, 2 or 3 in the ".ode" file or "ode.over" file. This is to make testing easier. I need to invent an example that would seem to require the higher number of sub-iterations.
June 3, 2012
Working on improving my HTML and omnisode documentation. I made documentation on the script files used to run omnisode.
June 2, 2012
Testing various numbers of sub-iterations with various increments:
1. H = 0.00001
2. H = 0.00005
3. H = 0.0001
4. H = 0.0002
5. H = 0.0005
6. H = 0.001
June 1, 2012
I am uncertain of something having to do with the solution of multiple (simultaneous) equations. Obviously some either have no solution or no unique solution. E.g.:
diff(y1,x,1) = diff(y2,x,1)
diff(y2,x,1) = diff(y1,x,1)
Has as solutions any differentiable equations where y1 and y2 differ by s constant.
When I first started I first tried sorting the equations - but this cannot always be done. Then I tried executing the code (atomall) in a loop, until the solution stopped changing much, but this was inadequate. Then I tried executing the loop the number of times there were equations, but I decided this was insufficient. I currently have the loop executing the sum of the orders of the equations times. I am not sure if this is adequate or not. Perhaps one needs to add the number of terms of the series to this. Of course, one still cannot arrive at a solution when there is none or a unique solution when there is none, but perhaps one can usually (except when the function is badly behaving for some other reason) obtain a solution in cases where there is, indeed, a unique solution.
I did some testing on this while watching the Celtics win.
May 31, 2012
Uploaded new test results with links to details.
Waiting till I have things better tested and documented before another release.
May 27, 2012
Testing a nonlinear equation.
Updated information on omnisode on Sourceforge and added latest tar file to its download page.
Msy 26, 2012
Ran automated test.
Also ran tests on systems of equations.
May 25, 2012
I have realized a slight problem with how I track the "revision" of the code in my table.html. Fairly often, especially when I reach a staple point, I run a script to save the three main pieces of ruby code. They are saved with numbers I specify, and the script prevents overwriting. It also saves the revision numbers in a file, used in the table.html when a test is recorded. Well, in testing I noticed when a test completed all the way to the end point (no timeout or maximum iterations reached), the cell containing the expected time remaining was not filled, and everything after it is shifted a cell to the left. So I worked fixing it. But I did not save my work after every attempt, so the "revision" codes only indicate the last time the code was saved - not the last time it was changed.
I've worked more on the html table The prior results are at html table 0001
May 24, 2012
I spent the last two days implementing automatic testing - constructing an html table automatically.
I haven't made any more changes to the algorithm. It will take a little time to decide how to package the automatic table generation.
May 23, 2012
The new tar file is omnisode0007.tar.gz which may be obtained at the Sourceforge download page for Omnisode.
There is one ode file format for both Maple and Maxima.
The ode file format was documented yesterday.
I briefly tested all ode files with maxima, and all seemed OK except complicated.ode and log.ode. Complicated seemed stuck. Maybe it was just very slow. Complicated.ode seemed OK with Maple. Also log.ode was not very accurate. But logs are difficult.
May 22, 2012
Maple: Sample generated Maple program
Maple: Sample Maple program results
Maxima: Sample generated Maxima program
Maple: Sample Maxima program results
I have written and tested making the ode files language independent. I have documented the ode files.
It will be tomorrow till I get the tar file made - less than a half hour.
May 21, 2012
More testing of new version (tar file from May 20).
See testrecords.html
I believe Maple is limited to about 32 digits precision. That's all I seem to get, and there is a negligible time difference in calculations assigning either 32 or 100 to "Digits".
May 20, 2012
After considerable effort trying to get my software to use maxima big floats, I realized that maxima does not support many functions (e.g. trig) for big floats. My software, thus, only supports 16 digits for maxima. I did, however, fix a few other things in trying. I have had time to do only limited testing, but I know of no problems, other than that the order of poles seems a little off, although the radius of convergence seems ok. Next I plan to do extensive testing and work on documentation.
omnisode0006.tar.gz was made.
May 19, 2012
It works for both Maple and Maxima now.
The biggest difficulty had to do with arrays in Maxima. It seems they are call by value.
I ended up using global arrays, and having the Ruby program generate code for each separate array.
I haven't had enough time to test enough, but the code is greatly improved!
The latest results are at:
testrecords.html
omnisode0004.tar.gz was made. I always worry that something will happen to me, before I finish.
May 15, 2012
It is working pretty well now.
I need to refamiliarize myself with sourceforge's tools but there is a tar file on my own web site for now:
omnitar0001.tar.gz was made. (removed as obsolete)
The testrecords.html has also been updated.
There is one sample output (and Maple echos the program) at Sample Solution
May 14, 2012
I have fixed a couple problems including multiple equations.
The results have been added to the testrecords.html
May 12, 2012
I am making some progress on this again the past few weeks.
I talk about this some on my blog Dennis' Daze
I am not releasing anything yet, but may start using the source version control tool soon.
A table of results of testing is at testrecords.html
December 12, 2010
I have abandoned this project.
There are problems I am unable to solve and my mind is not working well.
I intend to devote myself to philosophy - and may do some prolog programming for that.
These Ruby programs generate programs in Maple or Ruby to solve Systems of Ordinary Differential Equations. A long Taylor series method, pioneered by Prof. Y.F. Chang, who taught at the University of Nebraska in the late 1970's when I was a graduate student there, is used. The number of Taylor series terms can be specified in the problem file, though it is usually 30. The Maple version is gone as, at least for now, I do not have Maple available to me. I have started on a maxima version. In it I am taking time to remove some of the biggest inefficiencies. Also I may have found an error that I am trying to resolve. The Ruby version below would possibly also contain the error. But no one has reported any problems to me.
Coming soon here - a maxima version
There are problems with the Ruby and Maple versions which will have to wait until the maxima version is released.
Y. F. Chang's Draft Typescript(1978) on Taylor Series Chang on ATS
The Taylor series terms are used to calculate the values of the dependent variables, and also (optionally) the radius of convergence and order of any singularities. These can be used to reduce the size of the next increment. These can be used to reduce the size of the next increment. Increasing it in the middle results in large errors because each successive term is multiplied by increasing powers of the ratio of the new h to the old h. In the beginning I determine h by starting with a very small h and increase it until the estimated truncation error is as large as specified as permissible.
To visit the sourceforge project web page try Sode Project Page
For info on RubyApfp: RubyApfp Home Page

Systems of Ordinary Differential Equations

Table of Contents

Omnisode - What it is

Omnisode - Where it came from

Omnisode - Documentation

Omnisode - Latest Test Results

Omnisode - Known Bugs

Omnisode - Bugs Fixed

Omnisode - Download

Omnisode - Plans for the future

Omnisode - Items accomplished

Taylor Series - Bibliography

Notes