SODE
by Dennis J. Darland
djdarland@acm.org
dennis-darland73@augustana.net
SODE.icn is source to a Icon program to generate a Maple (or Icon) program to numerically solve either a single ordinary differential equation or a system of such equations. Whether "sodemaple.icn" or "sodeicon.icn" is included determines whether a Maple or Icon program is generated. sats.txt (for Maple) & sats.icn (for Icon) are libraries copied into the output files. The files in the archive may be placed in the directory of your choice, but sode must then be run in that directory & the problem files placed there as well.
The output maple code is found in "atomall.txt". (or "atomall.icn" for icon)
It is read into maple with
the command
read `c:\\path\\atomall.txt`;
Or compiled as usual wirh icon or unicon.
Then start the sode solver with "main();" within Maple.
It has been tested with Maple 7. There was a problem in Maple V 4.
The input file is found in a "---.ode"(in Maple) file (or "---.ide"(in Icon) file). There are many in this zip file. The format is below
Icon can be found at: http://www.cs.arizona.edu/icon/ Unicon can be found at: http://unicon.sourceforge.net/index.html Info on Maple is at: http://www.maplesoft.com/products/maple/
Format of .ode file:
ode spec
$
constants spec
$
initial values
$
exact_soln_y := proc(x) # analytic solutiion to be numerically compared. (or procedure exact_soln_y(x) for icon. ---
Any other functions( user defined)
$
The ode spec is of the form
diff(y,x,N) = x * sin(x);
The operations + - / * are supported.
The functions which may be used are
ln
exp
sin
cos
tan
sinh
cosh
tanh
arcsin
arccos
arctan
Si(not in Icon)
erf(not in Icon)
expt(a,b) (a to b power.)
Positive literals may also be used.
M1 can be used for -1.
diff (y,x,M) for M < N & M > 0
can be used for the lower derivatives of y.
SAMPLE .ide file
diff ( y1 , x , 6 ) = 3 * diff( y1, x , 1) - 2* diff( y2, x, 2); diff ( y2 , x , 3 ) = 2* diff( y1, x , 2) - 2 * diff( y3, x, 2); diff ( y3 , x , 3 ) = 2* diff( y1, x , 5) - 2 * y2; $
MAX_TERMS:=30; # Max Number of Taylor series terms MAX_TERMS:=11; # Min Number of Taylor series terms $ H := 0.0001; #initial H x_start := 3.0; #initial x x_end := 3.0003; #final x y1_init[1] := 1; #initial value of y1 y1_init[2] := 2; #initial value of y1' y1_init[3] := 0.7; #etc.
y1_init[4] := 1.7;
y1_init[5] := 0.6;
y1_init[6] := 0.9;
y2_init[1] := 9.6;
y2_init[2] := 6.5;
y2_init[3] := 8.5;
y3_init[1] := 6.8;
y3_init[2] := 6.5;
y3_init[3] := 3.5;
log10relerr := -6.0; #log base 10 of relative eror permittd default = -16; log10abserr := -6.0; #log base 10 of absolute eror permittd default = -16;
max_err := 10e-10; # IGNORED
HMIN := 0.00001; # minimum H
HMAX := 1.0; # maximum H adjust_h := 1; # whether to adjust H (1 is true. 0 is false.) look_poles := 1; # whether to estimate singularities(1 is true. 0 is false.)
# x_disp_incr := 0.1; # interval between which data is printed(Requires adjust_h) (Depreciated for disp_incr)
disp_incr :=0.1 # 0 for endpoints only -1 for "natural" H
MAX_ITER := 100; # Maximum number of iterations of main loop.
DELTA := 1.0e-10; # discrimination for display purposes.
SMALL_FLOAT := 1.0e-100; #default
SMALLISH_FLOAT := 1.0e-10; #default
LARGE_FLOAT := 1.0e100; #default
ALMOST_1 := 0.99; #default
Digits := 50; # Digits -- works only in maple.
$
procedure exact_soln_y1(x)
return 0.0; # this solution is unknown.
end
procedure exact_soln_y2(x)
return 0.0; # this solution is unknown.
end
procedure exact_soln_y3(x)
return 0.0; # this solution is unknown.
end
$