#BEGIN OUTFILE1
ALWAYS = 1
INFO = 2
DEBUGL = 3
DEBUGMASSIVE = 4
MAX_UNCHANGED = 10
# Begin Function number 3
def check_sign( x0 ,xf)
if (xf > x0) then # if number 1
ret = 1.0
else
ret = -1.0
end# end if 1
return ( ret)
end # End Function number 3
# Begin Function number 4
def est_size_answer()
min_size = $glob_large_float
if (omniabs($array_y2[1]) < min_size) then # if number 1
min_size = omniabs($array_y2[1])
omniout_float(ALWAYS,"min_size",32,min_size,32,"")
end# end if 1
if (omniabs($array_y1[1]) < min_size) then # if number 1
min_size = omniabs($array_y1[1])
omniout_float(ALWAYS,"min_size",32,min_size,32,"")
end# end if 1
if (min_size < 1.0) then # if number 1
min_size = 1.0
omniout_float(ALWAYS,"min_size",32,min_size,32,"")
end# end if 1
return ( min_size)
end # End Function number 4
# Begin Function number 5
def test_suggested_h()
max_value3 = 0.0
no_terms = $glob_max_terms
hn_div_ho = 0.5
hn_div_ho_2 = 0.25
hn_div_ho_3 = 0.125
omniout_float(ALWAYS,"hn_div_ho",32,hn_div_ho,32,"")
omniout_float(ALWAYS,"hn_div_ho_2",32,hn_div_ho_2,32,"")
omniout_float(ALWAYS,"hn_div_ho_3",32,hn_div_ho_3,32,"")
value3 = omniabs($array_y2[no_terms-3] + $array_y2[no_terms - 2] * hn_div_ho + $array_y2[no_terms - 1] * hn_div_ho_2 + $array_y2[no_terms] * hn_div_ho_3)
if (value3 > max_value3) then # if number 1
max_value3 = value3
omniout_float(ALWAYS,"value3",32,value3,32,"")
end# end if 1
value3 = omniabs($array_y1[no_terms-3] + $array_y1[no_terms - 2] * hn_div_ho + $array_y1[no_terms - 1] * hn_div_ho_2 + $array_y1[no_terms] * hn_div_ho_3)
if (value3 > max_value3) then # if number 1
max_value3 = value3
omniout_float(ALWAYS,"value3",32,value3,32,"")
end# end if 1
omniout_float(ALWAYS,"max_value3",32,max_value3,32,"")
return ( max_value3)
end # End Function number 5
# Begin Function number 6
def reached_interval()
if ($glob_check_sign * ($array_x[1]) >= $glob_check_sign * $glob_next_display) then # if number 1
ret = true
else
ret = false
end# end if 1
return(ret)
end # End Function number 6
# Begin Function number 7
def display_alot(iter)
#TOP DISPLAY ALOT
if (reached_interval()) then # if number 1
if (iter >= 0) then # if number 2
ind_var = $array_x[1]
omniout_float(ALWAYS,"x[1] ",33,ind_var,20," ")
analytic_val_y = exact_soln_y2(ind_var)
omniout_float(ALWAYS,"y2[1] (analytic) ",33,analytic_val_y,20," ")
term_no = 1
numeric_val = $array_y2[term_no]
abserr = omniabs(numeric_val - analytic_val_y)
omniout_float(ALWAYS,"y2[1] (numeric) ",33,numeric_val,20," ")
if (omniabs(analytic_val_y) != 0.0) then # if number 3
relerr = abserr*100.0/omniabs(analytic_val_y)
if (relerr > 0.0000000000000000000000000000000001) then # if number 4
$glob_good_digits = -trunc(log10(relerr)) + 2
else
$glob_good_digits = 16
end# end if 4
else
relerr = -1.0
$glob_good_digits = -1
end# end if 3
if ($glob_iter == 1) then # if number 3
$array_1st_rel_error[1] = relerr
else
$array_last_rel_error[1] = relerr
end# end if 3
omniout_float(ALWAYS,"absolute error ",4,abserr,20," ")
omniout_float(ALWAYS,"relative error ",4,relerr,20,"%")
omniout_int(INFO,"Correct digits ",32,$glob_good_digits,4," ")
omniout_float(ALWAYS,"h ",4,$glob_h,20," ")
analytic_val_y = exact_soln_y1(ind_var)
omniout_float(ALWAYS,"y1[1] (analytic) ",33,analytic_val_y,20," ")
term_no = 1
numeric_val = $array_y1[term_no]
abserr = omniabs(numeric_val - analytic_val_y)
omniout_float(ALWAYS,"y1[1] (numeric) ",33,numeric_val,20," ")
if (omniabs(analytic_val_y) != 0.0) then # if number 3
relerr = abserr*100.0/omniabs(analytic_val_y)
if (relerr > 0.0000000000000000000000000000000001) then # if number 4
$glob_good_digits = -trunc(log10(relerr)) + 2
else
$glob_good_digits = 16
end# end if 4
else
relerr = -1.0
$glob_good_digits = -1
end# end if 3
if ($glob_iter == 1) then # if number 3
$array_1st_rel_error[2] = relerr
else
$array_last_rel_error[2] = relerr
end# end if 3
omniout_float(ALWAYS,"absolute error ",4,abserr,20," ")
omniout_float(ALWAYS,"relative error ",4,relerr,20,"%")
omniout_int(INFO,"Correct digits ",32,$glob_good_digits,4," ")
omniout_float(ALWAYS,"h ",4,$glob_h,20," ")
end# end if 2
#BOTTOM DISPLAY ALOT
end# end if 1
end # End Function number 7
# Begin Function number 8
def adjust_for_pole(h_param)
#TOP ADJUST FOR POLE
hnew = h_param
$glob_normmax = $glob_small_float
if (omniabs($array_y2_higher[1][1]) > $glob_small_float) then # if number 1
tmp = omniabs($array_y2_higher[1][1])
if (tmp < $glob_normmax) then # if number 2
$glob_normmax = tmp
end# end if 2
end# end if 1
if (omniabs($array_y1_higher[1][1]) > $glob_small_float) then # if number 1
tmp = omniabs($array_y1_higher[1][1])
if (tmp < $glob_normmax) then # if number 2
$glob_normmax = tmp
end# end if 2
end# end if 1
if ($glob_look_poles and (omniabs($array_pole[1]) > $glob_small_float) and ($array_pole[1] != $glob_large_float)) then # if number 1
sz2 = $array_pole[1]/10.0
if (sz2 < hnew) then # if number 2
omniout_float(INFO,"$glob_h adjusted to ",20,h_param,12,"due to singularity.")
omniout_str(INFO,"Reached Optimal")
return(hnew)
end# end if 2
end# end if 1
if ( not $glob_reached_optimal_h) then # if number 1
$glob_reached_optimal_h = true
$glob_curr_iter_when_opt = $glob_current_iter
$glob_optimal_clock_start_sec = elapsed_time_seconds()
$glob_optimal_start = $array_x[1]
end# end if 1
hnew = sz2
#END block
return(hnew)
#BOTTOM ADJUST FOR POLE
end # End Function number 8
# Begin Function number 9
def prog_report(x_start,x_end)
#TOP PROGRESS REPORT
clock_sec1 = elapsed_time_seconds()
total_clock_sec = convfloat(clock_sec1) - convfloat($glob_orig_start_sec)
$glob_clock_sec = convfloat(clock_sec1) - convfloat($glob_clock_start_sec)
left_sec = convfloat($glob_max_sec) + convfloat($glob_orig_start_sec) - convfloat(clock_sec1)
expect_sec = comp_expect_sec(convfloat(x_end),convfloat(x_start),convfloat($array_x[1]) + convfloat($glob_h) ,convfloat( clock_sec1) - convfloat($glob_orig_start_sec))
opt_clock_sec = convfloat( clock_sec1) - convfloat($glob_optimal_clock_start_sec)
$glob_optimal_expect_sec = comp_expect_sec(convfloat(x_end),convfloat(x_start),convfloat($array_x[1]) +convfloat( $glob_h) ,convfloat( opt_clock_sec))
$glob_total_exp_sec = $glob_optimal_expect_sec + total_clock_sec
percent_done = comp_percent(convfloat(x_end),convfloat(x_start),convfloat($array_x[1]) + convfloat($glob_h))
$glob_percent_done = percent_done
omniout_str_noeol(INFO,"Total Elapsed Time ")
omniout_timestr(convfloat(total_clock_sec))
omniout_str_noeol(INFO,"Elapsed Time(since restart) ")
omniout_timestr(convfloat($glob_clock_sec))
if (convfloat(percent_done) < convfloat(100.0)) then # if number 1
omniout_str_noeol(INFO,"Expected Time Remaining ")
omniout_timestr(convfloat(expect_sec))
omniout_str_noeol(INFO,"Optimized Time Remaining ")
omniout_timestr(convfloat($glob_optimal_expect_sec))
omniout_str_noeol(INFO,"Expected Total Time ")
omniout_timestr(convfloat($glob_total_exp_sec))
end# end if 1
omniout_str_noeol(INFO,"Time to Timeout ")
omniout_timestr(convfloat(left_sec))
omniout_float(INFO, "Percent Done ",33,percent_done,4,"%")
#BOTTOM PROGRESS REPORT
end # End Function number 9
# Begin Function number 10
def check_for_pole()
#TOP CHECK FOR POLE
#IN RADII REAL EQ = 1
#Computes radius of convergence and r_order of pole from 3 adjacent Taylor series terms. EQUATUON NUMBER 1
#Applies to pole of arbitrary r_order on the real axis,
#Due to Prof. George Corliss.
n = $glob_max_terms
m = n - 5 - 1
while ((m >= 10) and ((omniabs($array_y2_higher[1][m]) < $glob_small_float * $glob_small_float) or (omniabs($array_y2_higher[1][m-1]) < $glob_small_float * $glob_small_float) or (omniabs($array_y2_higher[1][m-2]) < $glob_small_float * $glob_small_float ))) do # do number 2
m = m - 1
end# end do number 2
if (m > 10) then # if number 1
rm0 = $array_y2_higher[1][m]/$array_y2_higher[1][m-1]
rm1 = $array_y2_higher[1][m-1]/$array_y2_higher[1][m-2]
hdrc = convfloat(m)*rm0-convfloat(m-1)*rm1
if (omniabs(hdrc) > $glob_small_float * $glob_small_float) then # if number 2
rcs = $glob_h/hdrc
ord_no = (rm1*convfloat((m-2)*(m-2))-rm0*convfloat(m-3))/hdrc
$array_real_pole[1][1] = rcs
$array_real_pole[1][2] = ord_no
else
$array_real_pole[1][1] = $glob_large_float
$array_real_pole[1][2] = $glob_large_float
end# end if 2
else
$array_real_pole[1][1] = $glob_large_float
$array_real_pole[1][2] = $glob_large_float
end# end if 1
#BOTTOM RADII REAL EQ = 1
#IN RADII REAL EQ = 2
#Computes radius of convergence and r_order of pole from 3 adjacent Taylor series terms. EQUATUON NUMBER 2
#Applies to pole of arbitrary r_order on the real axis,
#Due to Prof. George Corliss.
n = $glob_max_terms
m = n - 1 - 1
while ((m >= 10) and ((omniabs($array_y1_higher[1][m]) < $glob_small_float * $glob_small_float) or (omniabs($array_y1_higher[1][m-1]) < $glob_small_float * $glob_small_float) or (omniabs($array_y1_higher[1][m-2]) < $glob_small_float * $glob_small_float ))) do # do number 2
m = m - 1
end# end do number 2
if (m > 10) then # if number 1
rm0 = $array_y1_higher[1][m]/$array_y1_higher[1][m-1]
rm1 = $array_y1_higher[1][m-1]/$array_y1_higher[1][m-2]
hdrc = convfloat(m)*rm0-convfloat(m-1)*rm1
if (omniabs(hdrc) > $glob_small_float * $glob_small_float) then # if number 2
rcs = $glob_h/hdrc
ord_no = (rm1*convfloat((m-2)*(m-2))-rm0*convfloat(m-3))/hdrc
$array_real_pole[2][1] = rcs
$array_real_pole[2][2] = ord_no
else
$array_real_pole[2][1] = $glob_large_float
$array_real_pole[2][2] = $glob_large_float
end# end if 2
else
$array_real_pole[2][1] = $glob_large_float
$array_real_pole[2][2] = $glob_large_float
end# end if 1
#BOTTOM RADII REAL EQ = 2
#TOP RADII COMPLEX EQ = 1
#Computes radius of convergence for complex conjugate pair of poles.
#from 6 adjacent Taylor series terms
#Also computes r_order of poles.
#Due to Manuel Prieto.
#With a correction by Dennis J. Darland
n = $glob_max_terms - 5 - 1
cnt = 0
while ((cnt < 5) and (n >= 10)) do # do number 2
if (omniabs($array_y2_higher[1][n]) > $glob_small_float) then # if number 1
cnt = cnt + 1
else
cnt = 0
end# end if 1
n = n - 1
end# end do number 2
m = n + cnt
if (m <= 10) then # if number 1
rad_c = $glob_large_float
ord_no = $glob_large_float
elsif
(((omniabs($array_y2_higher[1][m]) >= ($glob_large_float)) or (omniabs($array_y2_higher[1][m-1]) >=($glob_large_float)) or (omniabs($array_y2_higher[1][m-2]) >= ($glob_large_float)) or (omniabs($array_y2_higher[1][m-3]) >= ($glob_large_float)) or (omniabs($array_y2_higher[1][m-4]) >= ($glob_large_float)) or (omniabs($array_y2_higher[1][m-5]) >= ($glob_large_float))) or ((omniabs($array_y2_higher[1][m]) <= ($glob_small_float)) or (omniabs($array_y2_higher[1][m-1]) <=($glob_small_float)) or (omniabs($array_y2_higher[1][m-2]) <= ($glob_small_float)) or (omniabs($array_y2_higher[1][m-3]) <= ($glob_small_float)) or (omniabs($array_y2_higher[1][m-4]) <= ($glob_small_float)) or (omniabs($array_y2_higher[1][m-5]) <= ($glob_small_float)))) then # if number 2
rad_c = $glob_large_float
ord_no = $glob_large_float
else
rm0 = ($array_y2_higher[1][m])/($array_y2_higher[1][m-1])
rm1 = ($array_y2_higher[1][m-1])/($array_y2_higher[1][m-2])
rm2 = ($array_y2_higher[1][m-2])/($array_y2_higher[1][m-3])
rm3 = ($array_y2_higher[1][m-3])/($array_y2_higher[1][m-4])
rm4 = ($array_y2_higher[1][m-4])/($array_y2_higher[1][m-5])
nr1 = convfloat(m-1)*rm0 - 2.0*convfloat(m-2)*rm1 + convfloat(m-3)*rm2
nr2 = convfloat(m-2)*rm1 - 2.0*convfloat(m-3)*rm2 + convfloat(m-4)*rm3
dr1 = (-1.0)/rm1 + 2.0/rm2 - 1.0/rm3
dr2 = (-1.0)/rm2 + 2.0/rm3 - 1.0/rm4
ds1 = 3.0/rm1 - 8.0/rm2 + 5.0/rm3
ds2 = 3.0/rm2 - 8.0/rm3 + 5.0/rm4
if ((omniabs(nr1 * dr2 - nr2 * dr1) <= $glob_small_float) or (omniabs(dr1) <= $glob_small_float)) then # if number 3
rad_c = $glob_large_float
ord_no = $glob_large_float
else
if (omniabs(nr1*dr2 - nr2 * dr1) > $glob_small_float) then # if number 4
rcs = ((ds1*dr2 - ds2*dr1 +dr1*dr2)/(nr1*dr2 - nr2 * dr1))
#(Manuels) rcs = (ds1*dr2 - ds2*dr1)/(nr1*dr2 - nr2 * dr1)
ord_no = (rcs*nr1 - ds1)/(2.0*dr1) -convfloat(m)/2.0
if (omniabs(rcs) > $glob_small_float) then # if number 5
if (rcs > 0.0) then # if number 6
rad_c = sqrt(rcs) * omniabs($glob_h)
else
rad_c = $glob_large_float
end# end if 6
else
rad_c = $glob_large_float
ord_no = $glob_large_float
end# end if 5
else
rad_c = $glob_large_float
ord_no = $glob_large_float
end# end if 4
end# end if 3
$array_complex_pole[1][1] = rad_c
$array_complex_pole[1][2] = ord_no
end# end if 2
#BOTTOM RADII COMPLEX EQ = 1
#TOP RADII COMPLEX EQ = 2
#Computes radius of convergence for complex conjugate pair of poles.
#from 6 adjacent Taylor series terms
#Also computes r_order of poles.
#Due to Manuel Prieto.
#With a correction by Dennis J. Darland
n = $glob_max_terms - 1 - 1
cnt = 0
while ((cnt < 5) and (n >= 10)) do # do number 2
if (omniabs($array_y1_higher[1][n]) > $glob_small_float) then # if number 2
cnt = cnt + 1
else
cnt = 0
end# end if 2
n = n - 1
end# end do number 2
m = n + cnt
if (m <= 10) then # if number 2
rad_c = $glob_large_float
ord_no = $glob_large_float
elsif
(((omniabs($array_y1_higher[1][m]) >= ($glob_large_float)) or (omniabs($array_y1_higher[1][m-1]) >=($glob_large_float)) or (omniabs($array_y1_higher[1][m-2]) >= ($glob_large_float)) or (omniabs($array_y1_higher[1][m-3]) >= ($glob_large_float)) or (omniabs($array_y1_higher[1][m-4]) >= ($glob_large_float)) or (omniabs($array_y1_higher[1][m-5]) >= ($glob_large_float))) or ((omniabs($array_y1_higher[1][m]) <= ($glob_small_float)) or (omniabs($array_y1_higher[1][m-1]) <=($glob_small_float)) or (omniabs($array_y1_higher[1][m-2]) <= ($glob_small_float)) or (omniabs($array_y1_higher[1][m-3]) <= ($glob_small_float)) or (omniabs($array_y1_higher[1][m-4]) <= ($glob_small_float)) or (omniabs($array_y1_higher[1][m-5]) <= ($glob_small_float)))) then # if number 3
rad_c = $glob_large_float
ord_no = $glob_large_float
else
rm0 = ($array_y1_higher[1][m])/($array_y1_higher[1][m-1])
rm1 = ($array_y1_higher[1][m-1])/($array_y1_higher[1][m-2])
rm2 = ($array_y1_higher[1][m-2])/($array_y1_higher[1][m-3])
rm3 = ($array_y1_higher[1][m-3])/($array_y1_higher[1][m-4])
rm4 = ($array_y1_higher[1][m-4])/($array_y1_higher[1][m-5])
nr1 = convfloat(m-1)*rm0 - 2.0*convfloat(m-2)*rm1 + convfloat(m-3)*rm2
nr2 = convfloat(m-2)*rm1 - 2.0*convfloat(m-3)*rm2 + convfloat(m-4)*rm3
dr1 = (-1.0)/rm1 + 2.0/rm2 - 1.0/rm3
dr2 = (-1.0)/rm2 + 2.0/rm3 - 1.0/rm4
ds1 = 3.0/rm1 - 8.0/rm2 + 5.0/rm3
ds2 = 3.0/rm2 - 8.0/rm3 + 5.0/rm4
if ((omniabs(nr1 * dr2 - nr2 * dr1) <= $glob_small_float) or (omniabs(dr1) <= $glob_small_float)) then # if number 4
rad_c = $glob_large_float
ord_no = $glob_large_float
else
if (omniabs(nr1*dr2 - nr2 * dr1) > $glob_small_float) then # if number 5
rcs = ((ds1*dr2 - ds2*dr1 +dr1*dr2)/(nr1*dr2 - nr2 * dr1))
#(Manuels) rcs = (ds1*dr2 - ds2*dr1)/(nr1*dr2 - nr2 * dr1)
ord_no = (rcs*nr1 - ds1)/(2.0*dr1) -convfloat(m)/2.0
if (omniabs(rcs) > $glob_small_float) then # if number 6
if (rcs > 0.0) then # if number 7
rad_c = sqrt(rcs) * omniabs($glob_h)
else
rad_c = $glob_large_float
end# end if 7
else
rad_c = $glob_large_float
ord_no = $glob_large_float
end# end if 6
else
rad_c = $glob_large_float
ord_no = $glob_large_float
end# end if 5
end# end if 4
$array_complex_pole[2][1] = rad_c
$array_complex_pole[2][2] = ord_no
end# end if 3
#BOTTOM RADII COMPLEX EQ = 2
found_sing = 0
#TOP WHICH RADII EQ = 1
if (1 != found_sing and (($array_real_pole[1][1] == $glob_large_float) or ($array_real_pole[1][2] == $glob_large_float)) and (($array_complex_pole[1][1] != $glob_large_float) and ($array_complex_pole[1][2] != $glob_large_float)) and (($array_complex_pole[1][1] > 0.0) and ($array_complex_pole[1][2] > 0.0))) then # if number 3
$array_poles[1][1] = $array_complex_pole[1][1]
$array_poles[1][2] = $array_complex_pole[1][2]
found_sing = 1
$array_type_pole[1] = 2
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Complex estimate of poles used for equation 1")
end# end if 5
end# end if 4
end# end if 3
if (1 != found_sing and (($array_real_pole[1][1] != $glob_large_float) and ($array_real_pole[1][2] != $glob_large_float) and ($array_real_pole[1][1] > 0.0) and ($array_real_pole[1][2] > -1.0 * $glob_smallish_float) and (($array_complex_pole[1][1] == $glob_large_float) or ($array_complex_pole[1][2] == $glob_large_float) or ($array_complex_pole[1][1] <= 0.0 ) or ($array_complex_pole[1][2] <= 0.0)))) then # if number 3
$array_poles[1][1] = $array_real_pole[1][1]
$array_poles[1][2] = $array_real_pole[1][2]
found_sing = 1
$array_type_pole[1] = 1
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Real estimate of pole used for equation 1")
end# end if 5
end# end if 4
end# end if 3
if (1 != found_sing and ((($array_real_pole[1][1] == $glob_large_float) or ($array_real_pole[1][2] == $glob_large_float)) and (($array_complex_pole[1][1] == $glob_large_float) or ($array_complex_pole[1][2] == $glob_large_float)))) then # if number 3
$array_poles[1][1] = $glob_large_float
$array_poles[1][2] = $glob_large_float
found_sing = 1
$array_type_pole[1] = 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"NO POLE for equation 1")
end# end if 4
end# end if 3
if (1 != found_sing and (($array_real_pole[1][1] < $array_complex_pole[1][1]) and ($array_real_pole[1][1] > 0.0) and ($array_real_pole[1][2] > -1.0 * $glob_smallish_float))) then # if number 3
$array_poles[1][1] = $array_real_pole[1][1]
$array_poles[1][2] = $array_real_pole[1][2]
found_sing = 1
$array_type_pole[1] = 1
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Real estimate of pole used for equation 1")
end# end if 5
end# end if 4
end# end if 3
if (1 != found_sing and (($array_complex_pole[1][1] != $glob_large_float) and ($array_complex_pole[1][2] != $glob_large_float) and ($array_complex_pole[1][1] > 0.0) and ($array_complex_pole[1][2] > 0.0))) then # if number 3
$array_poles[1][1] = $array_complex_pole[1][1]
$array_poles[1][2] = $array_complex_pole[1][2]
$array_type_pole[1] = 2
found_sing = 1
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Complex estimate of poles used for equation 1")
end# end if 5
end# end if 4
end# end if 3
if (1 != found_sing ) then # if number 3
$array_poles[1][1] = $glob_large_float
$array_poles[1][2] = $glob_large_float
$array_type_pole[1] = 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"NO POLE for equation 1")
end# end if 4
end# end if 3
#BOTTOM WHICH RADII EQ = 1
#TOP WHICH RADII EQ = 2
if (2 != found_sing and (($array_real_pole[2][1] == $glob_large_float) or ($array_real_pole[2][2] == $glob_large_float)) and (($array_complex_pole[2][1] != $glob_large_float) and ($array_complex_pole[2][2] != $glob_large_float)) and (($array_complex_pole[2][1] > 0.0) and ($array_complex_pole[2][2] > 0.0))) then # if number 3
$array_poles[2][1] = $array_complex_pole[2][1]
$array_poles[2][2] = $array_complex_pole[2][2]
found_sing = 2
$array_type_pole[2] = 2
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Complex estimate of poles used for equation 2")
end# end if 5
end# end if 4
end# end if 3
if (2 != found_sing and (($array_real_pole[2][1] != $glob_large_float) and ($array_real_pole[2][2] != $glob_large_float) and ($array_real_pole[2][1] > 0.0) and ($array_real_pole[2][2] > -1.0 * $glob_smallish_float) and (($array_complex_pole[2][1] == $glob_large_float) or ($array_complex_pole[2][2] == $glob_large_float) or ($array_complex_pole[2][1] <= 0.0 ) or ($array_complex_pole[2][2] <= 0.0)))) then # if number 3
$array_poles[2][1] = $array_real_pole[2][1]
$array_poles[2][2] = $array_real_pole[2][2]
found_sing = 2
$array_type_pole[2] = 1
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Real estimate of pole used for equation 2")
end# end if 5
end# end if 4
end# end if 3
if (2 != found_sing and ((($array_real_pole[2][1] == $glob_large_float) or ($array_real_pole[2][2] == $glob_large_float)) and (($array_complex_pole[2][1] == $glob_large_float) or ($array_complex_pole[2][2] == $glob_large_float)))) then # if number 3
$array_poles[2][1] = $glob_large_float
$array_poles[2][2] = $glob_large_float
found_sing = 2
$array_type_pole[2] = 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"NO POLE for equation 2")
end# end if 4
end# end if 3
if (2 != found_sing and (($array_real_pole[2][1] < $array_complex_pole[2][1]) and ($array_real_pole[2][1] > 0.0) and ($array_real_pole[2][2] > -1.0 * $glob_smallish_float))) then # if number 3
$array_poles[2][1] = $array_real_pole[2][1]
$array_poles[2][2] = $array_real_pole[2][2]
found_sing = 2
$array_type_pole[2] = 1
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Real estimate of pole used for equation 2")
end# end if 5
end# end if 4
end# end if 3
if (2 != found_sing and (($array_complex_pole[2][1] != $glob_large_float) and ($array_complex_pole[2][2] != $glob_large_float) and ($array_complex_pole[2][1] > 0.0) and ($array_complex_pole[2][2] > 0.0))) then # if number 3
$array_poles[2][1] = $array_complex_pole[2][1]
$array_poles[2][2] = $array_complex_pole[2][2]
$array_type_pole[2] = 2
found_sing = 2
if ($glob_display_flag) then # if number 4
if (reached_interval()) then # if number 5
omniout_str(ALWAYS,"Complex estimate of poles used for equation 2")
end# end if 5
end# end if 4
end# end if 3
if (2 != found_sing ) then # if number 3
$array_poles[2][1] = $glob_large_float
$array_poles[2][2] = $glob_large_float
$array_type_pole[2] = 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"NO POLE for equation 2")
end# end if 4
end# end if 3
#BOTTOM WHICH RADII EQ = 2
$array_pole[1] = $glob_large_float
$array_pole[2] = $glob_large_float
#TOP WHICH RADIUS EQ = 1
if ($array_pole[1] > $array_poles[1][1]) then # if number 3
$array_pole[1] = $array_poles[1][1]
$array_pole[2] = $array_poles[1][2]
end# end if 3
#BOTTOM WHICH RADIUS EQ = 1
#TOP WHICH RADIUS EQ = 2
if ($array_pole[1] > $array_poles[2][1]) then # if number 3
$array_pole[1] = $array_poles[2][1]
$array_pole[2] = $array_poles[2][2]
end# end if 3
#BOTTOM WHICH RADIUS EQ = 2
#START ADJUST ALL SERIES
if ($array_pole[1] * $glob_ratio_of_radius < omniabs($glob_h)) then # if number 3
h_new = $array_pole[1] * $glob_ratio_of_radius
term = 1
ratio = 1.0
while (term <= $glob_max_terms) do # do number 2
$array_y2[term] = $array_y2[term]* ratio
$array_y2_higher[1][term] = $array_y2_higher[1][term]* ratio
$array_x[term] = $array_x[term]* ratio
$array_y1[term] = $array_y1[term]* ratio
$array_y1_higher[1][term] = $array_y1_higher[1][term]* ratio
$array_x[term] = $array_x[term]* ratio
ratio = ratio * h_new / omniabs($glob_h)
term = term + 1
end# end do number 2
$glob_h = h_new
end# end if 3
#BOTTOM ADJUST ALL SERIES
if (reached_interval()) then # if number 3
display_pole()
end# end if 3
end # End Function number 10
# Begin Function number 11
def get_norms()
if ( not $glob_initial_pass) then # if number 3
iii = 1
while (iii <= $glob_max_terms) do # do number 2
$array_norms[iii] = 0.0
iii = iii + 1
end# end do number 2
#TOP GET NORMS
iii = 1
while (iii <= $glob_max_terms) do # do number 2
if (omniabs($array_y2[iii]) > $array_norms[iii]) then # if number 4
$array_norms[iii] = omniabs($array_y2[iii])
end# end if 4
iii = iii + 1
end# end do number 2
iii = 1
while (iii <= $glob_max_terms) do # do number 2
if (omniabs($array_y1[iii]) > $array_norms[iii]) then # if number 4
$array_norms[iii] = omniabs($array_y1[iii])
end# end if 4
iii = iii + 1
end# end do number 2
#BOTTOM GET NORMS
end# end if 3
end # End Function number 11
# Begin Function number 12
def atomall()
#TOP ATOMALL
#END OUTFILE1
#BEGIN ATOMHDR1
#emit pre add CONST FULL $eq_no = 1 i = 1
$array_tmp1[1] = $array_const_0D0[1] + $array_y1[1]
#emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5
if ( not $array_y2_set_initial[1][6]) then # if number 1
if (1 <= $glob_max_terms) then # if number 2
temporary = $array_tmp1[1] * expt($glob_h , (5)) * factorial_3(0,5)
$array_y2[6] = temporary
$array_y2_higher[1][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y2_higher[2][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y2_higher[3][4] = temporary
temporary = temporary / $glob_h * (3.0)
$array_y2_higher[4][3] = temporary
temporary = temporary / $glob_h * (2.0)
$array_y2_higher[5][2] = temporary
temporary = temporary / $glob_h * (1.0)
$array_y2_higher[6][1] = temporary
end# end if 2
end# end if 1
kkk = 2
# emit pre mult FULL FULL $eq_no = 2 i = 1
$array_tmp3[1] = ($array_m1[1] * ($array_y2[1]))
#emit pre assign xxx $eq_no = 2 i = 1 $min_hdrs = 5
if ( not $array_y1_set_initial[2][2]) then # if number 1
if (1 <= $glob_max_terms) then # if number 2
temporary = $array_tmp3[1] * expt($glob_h , (1)) * factorial_3(0,1)
$array_y1[2] = temporary
$array_y1_higher[1][2] = temporary
temporary = temporary / $glob_h * (1.0)
$array_y1_higher[2][1] = temporary
end# end if 2
end# end if 1
kkk = 2
#END ATOMHDR1
#BEGIN ATOMHDR2
#emit pre add CONST FULL $eq_no = 1 i = 2
$array_tmp1[2] = $array_y1[2]
#emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5
if ( not $array_y2_set_initial[1][7]) then # if number 1
if (2 <= $glob_max_terms) then # if number 2
temporary = $array_tmp1[2] * expt($glob_h , (5)) * factorial_3(1,6)
$array_y2[7] = temporary
$array_y2_higher[1][7] = temporary
temporary = temporary / $glob_h * (6.0)
$array_y2_higher[2][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y2_higher[3][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y2_higher[4][4] = temporary
temporary = temporary / $glob_h * (3.0)
$array_y2_higher[5][3] = temporary
temporary = temporary / $glob_h * (2.0)
$array_y2_higher[6][2] = temporary
end# end if 2
end# end if 1
kkk = 3
# emit pre mult FULL FULL $eq_no = 2 i = 2
$array_tmp3[2] = ats(2,$array_m1,$array_y2,1)
#emit pre assign xxx $eq_no = 2 i = 2 $min_hdrs = 5
if ( not $array_y1_set_initial[2][3]) then # if number 1
if (2 <= $glob_max_terms) then # if number 2
temporary = $array_tmp3[2] * expt($glob_h , (1)) * factorial_3(1,2)
$array_y1[3] = temporary
$array_y1_higher[1][3] = temporary
temporary = temporary / $glob_h * (2.0)
$array_y1_higher[2][2] = temporary
end# end if 2
end# end if 1
kkk = 3
#END ATOMHDR2
#BEGIN ATOMHDR3
#emit pre add CONST FULL $eq_no = 1 i = 3
$array_tmp1[3] = $array_y1[3]
#emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5
if ( not $array_y2_set_initial[1][8]) then # if number 1
if (3 <= $glob_max_terms) then # if number 2
temporary = $array_tmp1[3] * expt($glob_h , (5)) * factorial_3(2,7)
$array_y2[8] = temporary
$array_y2_higher[1][8] = temporary
temporary = temporary / $glob_h * (7.0)
$array_y2_higher[2][7] = temporary
temporary = temporary / $glob_h * (6.0)
$array_y2_higher[3][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y2_higher[4][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y2_higher[5][4] = temporary
temporary = temporary / $glob_h * (3.0)
$array_y2_higher[6][3] = temporary
end# end if 2
end# end if 1
kkk = 4
# emit pre mult FULL FULL $eq_no = 2 i = 3
$array_tmp3[3] = ats(3,$array_m1,$array_y2,1)
#emit pre assign xxx $eq_no = 2 i = 3 $min_hdrs = 5
if ( not $array_y1_set_initial[2][4]) then # if number 1
if (3 <= $glob_max_terms) then # if number 2
temporary = $array_tmp3[3] * expt($glob_h , (1)) * factorial_3(2,3)
$array_y1[4] = temporary
$array_y1_higher[1][4] = temporary
temporary = temporary / $glob_h * (3.0)
$array_y1_higher[2][3] = temporary
end# end if 2
end# end if 1
kkk = 4
#END ATOMHDR3
#BEGIN ATOMHDR4
#emit pre add CONST FULL $eq_no = 1 i = 4
$array_tmp1[4] = $array_y1[4]
#emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5
if ( not $array_y2_set_initial[1][9]) then # if number 1
if (4 <= $glob_max_terms) then # if number 2
temporary = $array_tmp1[4] * expt($glob_h , (5)) * factorial_3(3,8)
$array_y2[9] = temporary
$array_y2_higher[1][9] = temporary
temporary = temporary / $glob_h * (8.0)
$array_y2_higher[2][8] = temporary
temporary = temporary / $glob_h * (7.0)
$array_y2_higher[3][7] = temporary
temporary = temporary / $glob_h * (6.0)
$array_y2_higher[4][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y2_higher[5][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y2_higher[6][4] = temporary
end# end if 2
end# end if 1
kkk = 5
# emit pre mult FULL FULL $eq_no = 2 i = 4
$array_tmp3[4] = ats(4,$array_m1,$array_y2,1)
#emit pre assign xxx $eq_no = 2 i = 4 $min_hdrs = 5
if ( not $array_y1_set_initial[2][5]) then # if number 1
if (4 <= $glob_max_terms) then # if number 2
temporary = $array_tmp3[4] * expt($glob_h , (1)) * factorial_3(3,4)
$array_y1[5] = temporary
$array_y1_higher[1][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y1_higher[2][4] = temporary
end# end if 2
end# end if 1
kkk = 5
#END ATOMHDR4
#BEGIN ATOMHDR5
#emit pre add CONST FULL $eq_no = 1 i = 5
$array_tmp1[5] = $array_y1[5]
#emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5
if ( not $array_y2_set_initial[1][10]) then # if number 1
if (5 <= $glob_max_terms) then # if number 2
temporary = $array_tmp1[5] * expt($glob_h , (5)) * factorial_3(4,9)
$array_y2[10] = temporary
$array_y2_higher[1][10] = temporary
temporary = temporary / $glob_h * (9.0)
$array_y2_higher[2][9] = temporary
temporary = temporary / $glob_h * (8.0)
$array_y2_higher[3][8] = temporary
temporary = temporary / $glob_h * (7.0)
$array_y2_higher[4][7] = temporary
temporary = temporary / $glob_h * (6.0)
$array_y2_higher[5][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y2_higher[6][5] = temporary
end# end if 2
end# end if 1
kkk = 6
# emit pre mult FULL FULL $eq_no = 2 i = 5
$array_tmp3[5] = ats(5,$array_m1,$array_y2,1)
#emit pre assign xxx $eq_no = 2 i = 5 $min_hdrs = 5
if ( not $array_y1_set_initial[2][6]) then # if number 1
if (5 <= $glob_max_terms) then # if number 2
temporary = $array_tmp3[5] * expt($glob_h , (1)) * factorial_3(4,5)
$array_y1[6] = temporary
$array_y1_higher[1][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y1_higher[2][5] = temporary
end# end if 2
end# end if 1
kkk = 6
#END ATOMHDR5
#BEGIN OUTFILE3
#Top Atomall While Loop-- outfile3
while (kkk <= $glob_max_terms) do # do number 1
#END OUTFILE3
#BEGIN OUTFILE4
#emit NOT FULL - FULL add $eq_no = 1
$array_tmp1[kkk] = $array_y1[kkk]
#emit assign $eq_no = 1
order_d = 5
if (kkk + order_d + 1 <= $glob_max_terms) then # if number 1
if ( not $array_y2_set_initial[1][kkk + order_d]) then # if number 2
temporary = $array_tmp1[kkk] * expt($glob_h , (order_d)) * factorial_3((kkk - 1),(kkk + order_d - 1))
$array_y2[kkk + order_d] = temporary
$array_y2_higher[1][kkk + order_d] = temporary
term = kkk + order_d - 1
adj2 = kkk + order_d - 1
adj3 = 2
while (term >= 1) do # do number 2
if (adj3 <= order_d + 1) then # if number 3
if (adj2 > 0) then # if number 4
temporary = temporary / $glob_h * convfp(adj2)
else
temporary = temporary
end# end if 4
$array_y2_higher[adj3][term] = temporary
end# end if 3
term = term - 1
adj2 = adj2 - 1
adj3 = adj3 + 1
end# end do number 2
end# end if 2
end# end if 1
#emit mult FULL FULL $eq_no = 2
$array_tmp3[kkk] = ats(kkk,$array_m1,$array_y2,1)
#emit assign $eq_no = 2
order_d = 1
if (kkk + order_d + 1 <= $glob_max_terms) then # if number 1
if ( not $array_y1_set_initial[2][kkk + order_d]) then # if number 2
temporary = $array_tmp3[kkk] * expt($glob_h , (order_d)) * factorial_3((kkk - 1),(kkk + order_d - 1))
$array_y1[kkk + order_d] = temporary
$array_y1_higher[1][kkk + order_d] = temporary
term = kkk + order_d - 1
adj2 = kkk + order_d - 1
adj3 = 2
while (term >= 1) do # do number 2
if (adj3 <= order_d + 1) then # if number 3
if (adj2 > 0) then # if number 4
temporary = temporary / $glob_h * convfp(adj2)
else
temporary = temporary
end# end if 4
$array_y1_higher[adj3][term] = temporary
end# end if 3
term = term - 1
adj2 = adj2 - 1
adj3 = adj3 + 1
end# end do number 2
end# end if 2
end# end if 1
kkk = kkk + 1
end# end do number 1
#BOTTOM ATOMALL
#END OUTFILE4
#BEGIN OUTFILE5
#BOTTOM ATOMALL ???
end # End Function number 12
#BEGIN ATS LIBRARY BLOCK
# Begin Function number 2
def omniout_str(iolevel,str)
if ($glob_iolevel >= iolevel) then # if number 1
puts str
end# end if 1
end # End Function number 2
# Begin Function number 3
def omniout_str_noeol(iolevel,str)
if ($glob_iolevel >= iolevel) then # if number 1
printf("%s", str)
end# end if 1
end # End Function number 3
# Begin Function number 4
def omniout_labstr(iolevel,label,str)
if ($glob_iolevel >= iolevel) then # if number 1
puts label + str
end# end if 1
end # End Function number 4
# Begin Function number 5
def omniout_float(iolevel,prelabel,prelen,value,vallen,postlabel)
if ($glob_iolevel >= iolevel) then # if number 1
puts prelabel.ljust(30) + value.to_s + postlabel
end# end if 1
end # End Function number 5
# Begin Function number 6
def omniout_int(iolevel,prelabel,prelen,value,vallen,postlabel)
if ($glob_iolevel >= iolevel) then # if number 1
puts prelabel.ljust(32) + value.to_s + postlabel
end# end if 1
end # End Function number 6
# Begin Function number 7
def omniout_float_arr(iolevel,prelabel,elemnt,prelen,value,vallen,postlabel)
if ($glob_iolevel >= iolevel) then # if number 1
puts (prelabel.ljust(12) + " " + elemnt.to_s + " " + value.to_s + " " + postlabel)
end# end if 1
end # End Function number 7
# Begin Function number 8
def dump_series(iolevel,dump_label,series_name,arr_series,numb)
if ($glob_iolevel >= iolevel) then # if number 1
i = 1
while (i <= numb) do
puts dump_label + series_name + "[" + i.to_s + "]" + arr_series[i].to_s
i += 1
end
end
end # End Function number 8
# Begin Function number 9
def dump_series_2(iolevel,dump_label,series_name2,arr_series2,numb,subnum,arr_x)
if ($glob_iolevel >= iolevel) then # if number 2
sub = 1;
while (sub <= subnum) do
i = 1;
while (i <= numb) do
puts dump_label + series_name2 + "[" + sub.to_s + "," + i.to_s + "]" + arr_series2[sub,i].to_s
end# end do number 0
sub += 1;
end# end do number -1
end# end if 2
end # End Function number 9
# Begin Function number 10
def cs_info(iolevel,str)
if ($glob_iolevel >= iolevel) then # if number 2
puts "cs_info " + str + " $glob_correct_start_flag = " , $glob_correct_start_flag.to_s + "$glob_h := " + $glob_h + "$glob_reached_optimal_h := " + $glob_reached_optimal_h
end# end if 2
end # End Function number 10
# Begin Function number 11
def logitem_time(fd,secs_in)
fd.printf("
")
if (secs_in >= 0) then # if number 2
years_int = trunc(secs_in / $glob_sec_in_year)
sec_temp = (secs_in % $glob_sec_in_year)
days_int = trunc(sec_temp / $glob_sec_in_day)
sec_temp = (sec_temp % $glob_sec_in_day)
hours_int = trunc(sec_temp / $glob_sec_in_hour)
sec_temp = (sec_temp % $glob_sec_in_hour)
minutes_int = trunc(sec_temp / $glob_sec_in_minute)
sec_int = (sec_temp % $glob_sec_in_minute)
if (years_int > 0) then # if number 3
fd.printf(years_int.to_s + " Years " + days_int.to_s + " Days " + hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds")
elsif
(days_int > 0) then # if number 4
fd.printf(days_int.to_s + " Days " + hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds")
elsif
(hours_int > 0) then # if number 5
fd.printf(hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds")
elsif
(minutes_int > 0) then # if number 6
fd.printf(minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds")
else
fd.printf(sec_int.to_s + " Seconds")
end# end if 6
else
fd.printf(" Unknown")
end# end if 5
fd.printf(" | ")
end # End Function number 11
# Begin Function number 12
def omniout_timestr(secs_in)
if (secs_in >= 0) then # if number 5
years_int = trunc(secs_in / $glob_sec_in_year)
sec_temp = (secs_in % $glob_sec_in_year)
days_int = trunc(sec_temp / $glob_sec_in_day)
sec_temp = (sec_temp % $glob_sec_in_day)
hours_int = trunc(sec_temp / $glob_sec_in_hour)
sec_temp = (sec_temp % $glob_sec_in_hour)
minutes_int = trunc(sec_temp / $glob_sec_in_minute)
sec_int = (sec_temp % $glob_sec_in_minute)
if (years_int > 0) then # if number 6
puts years_int.to_s + " Years " + days_int.to_s + " Days " + hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds"
elsif
(days_int > 0) then # if number 7
puts days_int.to_s + " Days " + hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds"
elsif
(hours_int > 0) then # if number 8
puts hours_int.to_s + " Hours " + minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds"
elsif
(minutes_int > 0) then # if number 9
puts minutes_int.to_s + " Minutes " + sec_int.to_s + " Seconds"
else
puts sec_int.to_s + " Seconds"
end# end if 9
else
puts " Unknown"
end# end if 8
end # End Function number 12
# Begin Function number 13
def ats(mmm_ats,arr_a,arr_b,jjj_ats)
ret_ats = 0.0
if (jjj_ats <= mmm_ats) then # if number 8
ma_ats = mmm_ats + 1
iii_ats = jjj_ats
while (iii_ats <= mmm_ats) do # do number -1
lll_ats = ma_ats - iii_ats
ret_ats = ret_ats + arr_a[iii_ats]*arr_b[lll_ats]
iii_ats = iii_ats + 1
end# end do number -1
end# end if 8
return ( ret_ats)
end # End Function number 13
# Begin Function number 14
def att(mmm_att,arr_aa,arr_bb,jjj_att)
ret_att = 0.0
if (jjj_att <= mmm_att) then # if number 8
ma_att = mmm_att + 2
iii_att = jjj_att
while (iii_att <= mmm_att) do # do number -1
lll_att = ma_att - iii_att
al_att = (lll_att - 1)
if (lll_att <= $glob_max_terms) then # if number 9
ret_att = ret_att + arr_aa[iii_att]*arr_bb[lll_att]* convfp(al_att)
end# end if 9
iii_att = iii_att + 1
end# end do number -1
ret_att = ret_att / convfp(mmm_att)
end# end if 8
return ( ret_att)
end # End Function number 14
# Begin Function number 15
def display_pole_debug(typ,radius,order2)
if (typ == 1) then # if number 8
omniout_str(ALWAYS,"Real")
else
omniout_str(ALWAYS,"Complex")
end# end if 8
omniout_float(ALWAYS,"DBG Radius of convergence ",4, radius,4," ")
omniout_float(ALWAYS,"DBG Order of pole ",4, order2,4," ")
end # End Function number 15
# Begin Function number 16
def display_pole()
if (($array_pole[1] != $glob_large_float) and ($array_pole[1] > 0.0) and ($array_pole[2] != $glob_large_float) and ($array_pole[2]> 0.0) and $glob_display_flag) then # if number 8
omniout_float(ALWAYS,"Radius of convergence ",4, $array_pole[1],4," ")
omniout_float(ALWAYS,"Order of pole ",4, $array_pole[2],4," ")
end# end if 8
end # End Function number 16
# Begin Function number 17
def logditto(file)
file.printf("")
file.printf("ditto")
file.printf(" | ")
end # End Function number 17
# Begin Function number 18
def logitem_integer(file,n)
file.printf("")
file.printf("%d",n)
file.printf(" | ")
end # End Function number 18
# Begin Function number 19
def logitem_str(file,str)
file.printf("")
file.printf(str)
file.printf(" | ")
end # End Function number 19
# Begin Function number 20
def logitem_good_digits(file,rel_error)
file.printf("")
if (rel_error != -1.0) then # if number 8
if (rel_error > + 0.0000000000000000000000000000000001) then # if number 9
good_digits = 1-trunc(log10(rel_error))
file.printf("%d",good_digits)
else
good_digits = 16
file.printf("%d",good_digits)
end# end if 9
else
file.printf("Unknown")
end# end if 8
file.printf(" | ")
end # End Function number 20
# Begin Function number 21
def log_revs(file,revs)
file.printf(revs.to_s)
end # End Function number 21
# Begin Function number 22
def logitem_float(file,x)
file.printf("")
file.printf(x.to_s)
file.printf(" | ")
end # End Function number 22
# Begin Function number 23
def logitem_pole(file,pole)
file.printf("")
if pole == 0 then # if number 8
file.printf("NA")
elsif
pole == 1 then # if number 9
file.printf("Real")
elsif
pole == 2 then # if number 10
file.printf("Complex")
else
file.printf("No Pole")
end# end if 10
file.printf(" | ")
end # End Function number 23
# Begin Function number 24
def logstart(file)
file.printf("")
end # End Function number 24
# Begin Function number 25
def logend(file)
file.printf("
")
end # End Function number 25
# Begin Function number 26
def chk_data()
errflag = false
if (($glob_max_terms < 15) or ($glob_max_terms > 512)) then # if number 10
omniout_str(ALWAYS,"Illegal max_terms = -- Using 30")
$glob_max_terms = 30
end# end if 10
if ($glob_max_iter < 2) then # if number 10
omniout_str(ALWAYS,"Illegal max_iter")
errflag = true
end# end if 10
if (errflag) then # if number 10
end# end if 10
end # End Function number 26
# Begin Function number 27
def comp_expect_sec(t_end2,t_start2,t2,clock_sec2)
ms2 = clock_sec2
sub1 = (t_end2-t_start2)
sub2 = (t2-t_start2)
if (sub1 == 0.0) then # if number 10
sec_left = 0.0
else
if (sub2 > 0.0) then # if number 11
rrr = (sub1/sub2)
sec_left = rrr * ms2 - ms2
else
sec_left = 0.0
end# end if 11
end# end if 10
return ( sec_left)
end # End Function number 27
# Begin Function number 28
def comp_percent(t_end2,t_start2, t2)
sub1 = (t_end2-t_start2)
sub2 = (t2-t_start2)
if (sub2 > $glob_small_float) then # if number 10
rrr = (100.0*sub2)/sub1
else
rrr = 0.0
end# end if 10
return ( rrr)
end # End Function number 28
# Begin Function number 29
def factorial_2(nnn)
if (nnn > 0) then # if number 10
return ( nnn * factorial_2(nnn - 1))
else
return ( 1.0)
end# end if 10
end # End Function number 29
# Begin Function number 30
def factorial_1(nnn)
if (nnn <= $glob_max_terms) then # if number 10
if ($array_fact_1[nnn] == 0) then # if number 11
ret = factorial_2(nnn)
$array_fact_1[nnn] = ret
else
ret = $array_fact_1[nnn]
end# end if 11
else
ret = factorial_2(nnn)
end# end if 10
return ( ret)
end # End Function number 30
# Begin Function number 31
def factorial_3(mmm,nnn)
if ((nnn <= $glob_max_terms) and (mmm <= $glob_max_terms)) then # if number 10
if ($array_fact_2[mmm][nnn] == 0) then # if number 11
ret = factorial_1(mmm)/factorial_1(nnn)
$array_fact_2[mmm][nnn] = ret
else
ret = $array_fact_2[mmm][nnn]
end# end if 11
else
ret = factorial_2(mmm)/factorial_2(nnn)
end# end if 10
return ( ret)
end # End Function number 31
# Begin Function number 32
def convfp(mmm)
return(mmm.to_f)
end # End Function number 32
# Begin Function number 33
def convfloat(mmm)
return(mmm.to_f)
end # End Function number 33
# Begin Function number 34
def elapsed_time_seconds()
t = Time.now
return(t.to_i)
end # End Function number 34
# Begin Function number 35
def expt(x,y)
return(x**y)
end # End Function number 35
# Begin Function number 36
def Si(x)
return(0.0)
end # End Function number 36
# Begin Function number 37
def Ci(x)
return(0.0)
end # End Function number 37
# Begin Function number 38
def abs(x)
return(x.abs)
end # End Function number 38
# Begin Function number 39
def omniabs(x)
return(x.abs)
end # End Function number 39
# Begin Function number 40
def exp(x)
return(Math.exp(x))
end # End Function number 40
# Begin Function number 41
def trunc(x)
return(x.to_i)
end # End Function number 41
# Begin Function number 42
def floor(x)
return(x.floor)
end # End Function number 42
# Begin Function number 43
def sin(x)
return(Math.sin(x))
end # End Function number 43
# Begin Function number 44
def cos(x)
return(Math.cos(x))
end # End Function number 44
# Begin Function number 45
def tan(x)
return(Math.tan(x))
end # End Function number 45
# Begin Function number 46
def arccos(x)
return(Math.acos(x))
end # End Function number 46
# Begin Function number 47
def arccosh(x)
return(Math.acosh(x))
end # End Function number 47
# Begin Function number 48
def arcsin(x)
return(Math.asin(x))
end # End Function number 48
# Begin Function number 49
def arcsinh(x)
return(Math.asinh(x))
end # End Function number 49
# Begin Function number 50
def arctan(x)
return(Math.atan(x))
end # End Function number 50
# Begin Function number 51
def arctanh(x)
return(Math.atanh(x))
end # End Function number 51
# Begin Function number 52
def cosh(x)
return(Math.cosh(x))
end # End Function number 52
# Begin Function number 53
def erf(x)
return(Math.erf(x))
end # End Function number 53
# Begin Function number 54
def ln(x)
return(Math.log(x))
end # End Function number 54
# Begin Function number 55
def log10(x)
return(Math.log10(x))
end # End Function number 55
# Begin Function number 56
def sinh(x)
return(Math.sinh(x))
end # End Function number 56
# Begin Function number 57
def tanh(x)
return(Math.tanh(x))
end # End Function number 57
# Begin Function number 58
def sqrt(x)
return(Math.sqrt(x))
end # End Function number 58
# Begin Function number 59
def array2d(op3,op4)
i = 0
x1 = Array.new(op3.to_i + 1)
while i <= op3.to_i + 1 do # do number -1
x1[i] = Array.new(op4.to_i + 1)
i += 1
end# end do number -1
return x1
end # End Function number 59
# Begin Function number 60
def estimated_needed_step_error(x_start,x_end,estimated_h,estimated_answer)
omniout_float(ALWAYS,"$glob_desired_digits_correct",32,$glob_desired_digits_correct,32,"")
desired_abs_gbl_error = expt(10.0,- $glob_desired_digits_correct) * omniabs(estimated_answer)
omniout_float(ALWAYS,"desired_abs_gbl_error",32,desired_abs_gbl_error,32,"")
range = (x_end - x_start)
omniout_float(ALWAYS,"range",32,range,32,"")
estimated_steps = range / estimated_h
omniout_float(ALWAYS,"estimated_steps",32,estimated_steps,32,"")
step_error = omniabs(desired_abs_gbl_error / estimated_steps)
omniout_float(ALWAYS,"step_error",32,step_error,32,"")
return ( (step_error))
end # End Function number 60
#END ATS LIBRARY BLOCK
#BEGIN USER DEF BLOCK
#BEGIN USER DEF BLOCK
def exact_soln_y1 (x)
return( cos(x))
end
def exact_soln_y2 (x)
return( sin(x))
end
def exact_soln_y2p (x)
return( cos(x))
end
def exact_soln_y2pp (x)
return( -sin(x))
end
def exact_soln_y2ppp (x)
return( -cos(x))
end
def exact_soln_y2pppp (x)
return( sin(x))
end
#END USER DEF BLOCK
#END USER DEF BLOCK
#END OUTFILE5
# Begin Function number 2
def main()
#BEGIN OUTFIEMAIN
$glob_iolevel = INFO
$glob_check_sign = 1.0
$glob_desired_digits_correct = 8.0
$glob_max_value3 = 0.0
$glob_ratio_of_radius = 0.01
$glob_percent_done = 0.0
$glob_subiter_method = 3
$glob_total_exp_sec = 0.1
$glob_optimal_expect_sec = 0.1
$glob_html_log = true
$glob_good_digits = 0
$glob_max_opt_iter = 10
$glob_dump = false
$glob_djd_debug = true
$glob_display_flag = true
$glob_djd_debug2 = true
$glob_sec_in_minute = 60
$glob_min_in_hour = 60
$glob_hours_in_day = 24
$glob_days_in_year = 365
$glob_sec_in_hour = 3600
$glob_sec_in_day = 86400
$glob_sec_in_year = 31536000
$glob_almost_1 = 0.9990
$glob_clock_sec = 0.0
$glob_clock_start_sec = 0.0
$glob_not_yet_finished = true
$glob_initial_pass = true
$glob_not_yet_start_msg = true
$glob_reached_optimal_h = false
$glob_optimal_done = false
$glob_disp_incr = 0.1
$glob_h = 0.1
$glob_max_h = 0.1
$glob_large_float = 9.0e100
$glob_last_good_h = 0.1
$glob_look_poles = false
$glob_neg_h = false
$glob_display_interval = 0.0
$glob_next_display = 0.0
$glob_dump_analytic = false
$glob_abserr = 0.1e-10
$glob_relerr = 0.1e-10
$glob_max_hours = 0.0
$glob_max_iter = 1000
$glob_max_rel_trunc_err = 0.1e-10
$glob_max_trunc_err = 0.1e-10
$glob_no_eqs = 0
$glob_optimal_clock_start_sec = 0.0
$glob_optimal_start = 0.0
$glob_small_float = 0.1e-200
$glob_smallish_float = 0.1e-100
$glob_unchanged_h_cnt = 0
$glob_warned = false
$glob_warned2 = false
$glob_max_sec = 10000.0
$glob_orig_start_sec = 0.0
$glob_start = 0
$glob_curr_iter_when_opt = 0
$glob_current_iter = 0
$glob_iter = 0
$glob_normmax = 0.0
$glob_max_minutes = 0.0
#Write Set Defaults
$glob_orig_start_sec = elapsed_time_seconds()
$glob_curr_iter_when_opt = 0
$glob_display_flag = true
$glob_no_eqs = 2
$glob_iter = -1
opt_iter = -1
$glob_max_iter = 50000
$glob_max_hours = 0.0
$glob_max_minutes = 15.0
omniout_str(ALWAYS,"##############ECHO OF PROBLEM#################")
omniout_str(ALWAYS,"##############temp/mtest7postode.ode#################")
omniout_str(ALWAYS,"diff ( y2 , x , 5 ) = y1 ;")
omniout_str(ALWAYS,"diff ( y1 , x , 1 ) = m1 * y2 ;")
omniout_str(ALWAYS,"!")
omniout_str(ALWAYS,"#BEGIN FIRST INPUT BLOCK")
omniout_str(ALWAYS,"# Digits:=32; ELIMINATED in preodein.rb")
omniout_str(ALWAYS,"max_terms=30 ")
omniout_str(ALWAYS,"!")
omniout_str(ALWAYS,"#END FIRST INPUT BLOCK")
omniout_str(ALWAYS,"#BEGIN SECOND INPUT BLOCK")
omniout_str(ALWAYS,"x_start=0.0 ")
omniout_str(ALWAYS,"x_end=5.0 ")
omniout_str(ALWAYS,"$array_y1_init[0 + 1] = exact_soln_y1(x_start)")
omniout_str(ALWAYS,"$array_y2_init[0 + 1] = exact_soln_y2(x_start)")
omniout_str(ALWAYS,"$array_y2_init[1 + 1] = exact_soln_y2p(x_start)")
omniout_str(ALWAYS,"$array_y2_init[2 + 1] = exact_soln_y2pp(x_start)")
omniout_str(ALWAYS,"$array_y2_init[4 + 1] = exact_soln_y2pppp(x_start)")
omniout_str(ALWAYS,"$glob_look_poles=true ")
omniout_str(ALWAYS,"$glob_max_iter=20 ")
omniout_str(ALWAYS,"#END SECOND INPUT BLOCK")
omniout_str(ALWAYS,"#BEGIN OVERRIDE BLOCK")
omniout_str(ALWAYS,"$glob_desired_digits_correct=10 ")
omniout_str(ALWAYS,"$glob_display_interval=0.001 ")
omniout_str(ALWAYS,"$glob_look_poles=true ")
omniout_str(ALWAYS,"$glob_max_iter=10000000 ")
omniout_str(ALWAYS,"$glob_max_minutes=3 ")
omniout_str(ALWAYS,"$glob_subiter_method=3 ")
omniout_str(ALWAYS,"#END OVERRIDE BLOCK")
omniout_str(ALWAYS,"!")
omniout_str(ALWAYS,"#BEGIN USER DEF BLOCK")
omniout_str(ALWAYS,"def exact_soln_y1 (x)")
omniout_str(ALWAYS,"return( cos(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"def exact_soln_y2 (x)")
omniout_str(ALWAYS,"return( sin(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"def exact_soln_y2p (x)")
omniout_str(ALWAYS,"return( cos(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"def exact_soln_y2pp (x)")
omniout_str(ALWAYS,"return( -sin(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"def exact_soln_y2ppp (x)")
omniout_str(ALWAYS,"return( -cos(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"def exact_soln_y2pppp (x)")
omniout_str(ALWAYS,"return( sin(x)) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"")
omniout_str(ALWAYS,"#END USER DEF BLOCK")
omniout_str(ALWAYS,"#######END OF ECHO OF PROBLEM#################")
$glob_unchanged_h_cnt = 0
$glob_warned = false
$glob_warned2 = false
$glob_small_float = 1.0e-200
$glob_smallish_float = 1.0e-64
$glob_large_float = 1.0e100
$glob_almost_1 = 0.99
#BEGIN FIRST INPUT BLOCK
#BEGIN FIRST INPUT BLOCK
# Digits:=32; ELIMINATED in preodein.rb
max_terms=30
#END FIRST INPUT BLOCK
#START OF INITS AFTER INPUT BLOCK
$glob_max_terms = max_terms
$glob_html_log = true
#END OF INITS AFTER INPUT BLOCK
$array_y2_init= Array.new(max_terms + 1)
$array_y1_init= Array.new(max_terms + 1)
$array_norms= Array.new(max_terms + 1)
$array_fact_1= Array.new(max_terms + 1)
$array_pole= Array.new(max_terms + 1)
$array_1st_rel_error= Array.new(max_terms + 1)
$array_last_rel_error= Array.new(max_terms + 1)
$array_type_pole= Array.new(max_terms + 1)
$array_y2= Array.new(max_terms + 1)
$array_x= Array.new(max_terms + 1)
$array_y1= Array.new(max_terms + 1)
$array_tmp0= Array.new(max_terms + 1)
$array_tmp1= Array.new(max_terms + 1)
$array_tmp2= Array.new(max_terms + 1)
$array_tmp3= Array.new(max_terms + 1)
$array_m1= Array.new(max_terms + 1)
$array_y2_higher = array2d(6 + 1 ,max_terms + 1 )
$array_y2_higher_work = array2d(6 + 1 ,max_terms + 1 )
$array_y2_higher_work2 = array2d(6 + 1 ,max_terms + 1 )
$array_y2_set_initial = array2d(3 + 1 ,max_terms + 1 )
$array_y1_higher = array2d(2 + 1 ,max_terms + 1 )
$array_y1_higher_work = array2d(2 + 1 ,max_terms + 1 )
$array_y1_higher_work2 = array2d(2 + 1 ,max_terms + 1 )
$array_y1_set_initial = array2d(3 + 1 ,max_terms + 1 )
$array_poles = array2d(2 + 1 ,3 + 1 )
$array_real_pole = array2d(2 + 1 ,3 + 1 )
$array_complex_pole = array2d(2 + 1 ,3 + 1 )
$array_fact_2 = array2d(max_terms + 1 ,max_terms + 1 )
term = 1
while (term <= max_terms) do # do number 2
$array_y2_init[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_y1_init[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_norms[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_fact_1[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_pole[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_1st_rel_error[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_last_rel_error[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_type_pole[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_y2[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_x[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_y1[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp0[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp1[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp2[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp3[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_m1[term] = 0.0
term = term + 1
end# end do number 2
ord = 1
while (ord <=6) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y2_higher[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=6) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y2_higher_work[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=6) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y2_higher_work2[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=3) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y2_set_initial[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y1_higher[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y1_higher_work[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y1_higher_work2[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=3) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y1_set_initial[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= 3) do # do number 3
$array_poles[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= 3) do # do number 3
$array_real_pole[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=2) do # do number 2
term = 1
while (term <= 3) do # do number 3
$array_complex_pole[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
ord = 1
while (ord <=max_terms) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_fact_2[ord][term] = 0.0
term = term + 1
end# end do number 3
ord = ord + 1
end# end do number 2
#BEGIN ARRAYS DEFINED AND INITIALIZATED
$array_y2 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_y2[term] = 0.0
term = term + 1
end# end do number 2
$array_x = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_x[term] = 0.0
term = term + 1
end# end do number 2
$array_y1 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_y1[term] = 0.0
term = term + 1
end# end do number 2
$array_m1 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_m1[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp0 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp0[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp1 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp1[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp2 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp2[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp3 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp3[term] = 0.0
term = term + 1
end# end do number 2
$array_const_5 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_const_5[term] = 0.0
term = term + 1
end# end do number 2
$array_const_5[1] = 5
$array_const_0D0 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_const_0D0[term] = 0.0
term = term + 1
end# end do number 2
$array_const_0D0[1] = 0.0
$array_const_1 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_const_1[term] = 0.0
term = term + 1
end# end do number 2
$array_const_1[1] = 1
$array_m1 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms) do # do number 2
$array_m1[term] = 0.0
term = term + 1
end# end do number 2
$array_m1[1] = -1.0
#END ARRAYS DEFINED AND INITIALIZATED
#Initing Factorial Tables
iiif = 0
while (iiif <= $glob_max_terms) do # do number 2
jjjf = 0
while (jjjf <= $glob_max_terms) do # do number 3
$array_fact_1[iiif] = 0
$array_fact_2[iiif][jjjf] = 0
jjjf = jjjf + 1
end# end do number 3
iiif = iiif + 1
end# end do number 2
#Done Initing Factorial Tables
#TOP SECOND INPUT BLOCK
#BEGIN SECOND INPUT BLOCK
#END FIRST INPUT BLOCK
#BEGIN SECOND INPUT BLOCK
x_start=0.0
x_end=5.0
$array_y1_init[0 + 1] = exact_soln_y1(x_start)
$array_y2_init[0 + 1] = exact_soln_y2(x_start)
$array_y2_init[1 + 1] = exact_soln_y2p(x_start)
$array_y2_init[2 + 1] = exact_soln_y2pp(x_start)
$array_y2_init[4 + 1] = exact_soln_y2pppp(x_start)
$glob_look_poles=true
$glob_max_iter=20
#END SECOND INPUT BLOCK
#BEGIN OVERRIDE BLOCK
$glob_desired_digits_correct=10
$glob_display_interval=0.001
$glob_look_poles=true
$glob_max_iter=10000000
$glob_max_minutes=3
$glob_subiter_method=3
#END OVERRIDE BLOCK
#END SECOND INPUT BLOCK
#BEGIN INITS AFTER SECOND INPUT BLOCK
$glob_last_good_h = $glob_h
$glob_max_terms = max_terms
$glob_max_sec = convfloat(60.0) * convfloat($glob_max_minutes) + convfloat(3600.0) * convfloat($glob_max_hours)
if ($glob_h > 0.0) then # if number 1
$glob_neg_h = false
$glob_display_interval = omniabs($glob_display_interval)
else
$glob_neg_h = true
$glob_display_interval = -omniabs($glob_display_interval)
end# end if 1
chk_data()
#AFTER INITS AFTER SECOND INPUT BLOCK
$array_y2_set_initial[1][1] = true
$array_y2_set_initial[1][2] = true
$array_y2_set_initial[1][3] = true
$array_y2_set_initial[1][4] = false
$array_y2_set_initial[1][5] = true
$array_y2_set_initial[1][6] = false
$array_y2_set_initial[1][7] = false
$array_y2_set_initial[1][8] = false
$array_y2_set_initial[1][9] = false
$array_y2_set_initial[1][10] = false
$array_y2_set_initial[1][11] = false
$array_y2_set_initial[1][12] = false
$array_y2_set_initial[1][13] = false
$array_y2_set_initial[1][14] = false
$array_y2_set_initial[1][15] = false
$array_y2_set_initial[1][16] = false
$array_y2_set_initial[1][17] = false
$array_y2_set_initial[1][18] = false
$array_y2_set_initial[1][19] = false
$array_y2_set_initial[1][20] = false
$array_y2_set_initial[1][21] = false
$array_y2_set_initial[1][22] = false
$array_y2_set_initial[1][23] = false
$array_y2_set_initial[1][24] = false
$array_y2_set_initial[1][25] = false
$array_y2_set_initial[1][26] = false
$array_y2_set_initial[1][27] = false
$array_y2_set_initial[1][28] = false
$array_y2_set_initial[1][29] = false
$array_y2_set_initial[1][30] = false
$array_y1_set_initial[2][1] = true
$array_y1_set_initial[2][2] = false
$array_y1_set_initial[2][3] = false
$array_y1_set_initial[2][4] = false
$array_y1_set_initial[2][5] = false
$array_y1_set_initial[2][6] = false
$array_y1_set_initial[2][7] = false
$array_y1_set_initial[2][8] = false
$array_y1_set_initial[2][9] = false
$array_y1_set_initial[2][10] = false
$array_y1_set_initial[2][11] = false
$array_y1_set_initial[2][12] = false
$array_y1_set_initial[2][13] = false
$array_y1_set_initial[2][14] = false
$array_y1_set_initial[2][15] = false
$array_y1_set_initial[2][16] = false
$array_y1_set_initial[2][17] = false
$array_y1_set_initial[2][18] = false
$array_y1_set_initial[2][19] = false
$array_y1_set_initial[2][20] = false
$array_y1_set_initial[2][21] = false
$array_y1_set_initial[2][22] = false
$array_y1_set_initial[2][23] = false
$array_y1_set_initial[2][24] = false
$array_y1_set_initial[2][25] = false
$array_y1_set_initial[2][26] = false
$array_y1_set_initial[2][27] = false
$array_y1_set_initial[2][28] = false
$array_y1_set_initial[2][29] = false
$array_y1_set_initial[2][30] = false
#BEGIN OPTIMIZE CODE
omniout_str(ALWAYS,"START of Optimize")
#Start Series -- INITIALIZE FOR OPTIMIZE
$glob_check_sign = check_sign(x_start,x_end)
$glob_h = check_sign(x_start,x_end)
if ($glob_display_interval < $glob_h) then # if number 3
$glob_h = $glob_display_interval
end# end if 3
if ($glob_max_h < $glob_h) then # if number 3
$glob_h = $glob_max_h
end# end if 3
found_h = -1.0
best_h = 0.0
min_value = $glob_large_float
est_answer = est_size_answer()
opt_iter = 1
while ((opt_iter <= 20) and (found_h < 0.0)) do # do number 2
omniout_int(ALWAYS,"opt_iter",32,opt_iter,4,"")
$array_x[1] = x_start
$array_x[2] = $glob_h
$glob_next_display = x_start
order_diff = 5
#Start Series $array_y2
term_no = 1
while (term_no <= order_diff) do # do number 3
$array_y2[term_no] = $array_y2_init[term_no] * expt($glob_h , (term_no - 1)) / factorial_1(term_no - 1)
term_no = term_no + 1
end# end do number 3
rows = order_diff
r_order = 1
while (r_order <= rows) do # do number 3
term_no = 1
while (term_no <= (rows - r_order + 1)) do # do number 4
it = term_no + r_order - 1
$array_y2_higher[r_order][term_no] = $array_y2_init[it]* expt($glob_h , (term_no - 1)) / ((factorial_1(term_no - 1)))
term_no = term_no + 1
end# end do number 4
r_order = r_order + 1
end# end do number 3
order_diff = 1
#Start Series $array_y1
term_no = 1
while (term_no <= order_diff) do # do number 3
$array_y1[term_no] = $array_y1_init[term_no] * expt($glob_h , (term_no - 1)) / factorial_1(term_no - 1)
term_no = term_no + 1
end# end do number 3
rows = order_diff
r_order = 1
while (r_order <= rows) do # do number 3
term_no = 1
while (term_no <= (rows - r_order + 1)) do # do number 4
it = term_no + r_order - 1
$array_y1_higher[r_order][term_no] = $array_y1_init[it]* expt($glob_h , (term_no - 1)) / ((factorial_1(term_no - 1)))
term_no = term_no + 1
end# end do number 4
r_order = r_order + 1
end# end do number 3
if ($glob_subiter_method == 1 ) then # if number 3
atomall()
elsif
($glob_subiter_method == 2 ) then # if number 4
subiter = 1
while (subiter <= 6) do # do number 3
atomall()
subiter = subiter + 1
end# end do number 3
else
subiter = 1
while (subiter <= 6 + $glob_max_terms) do # do number 3
atomall()
subiter = subiter + 1
end# end do number 3
end# end if 4
est_needed_step_err = estimated_needed_step_error(x_start,x_end,$glob_h,est_answer)
omniout_float(ALWAYS,"est_needed_step_err",32,est_needed_step_err,16,"")
value3 = test_suggested_h()
omniout_float(ALWAYS,"value3",32,value3,32,"")
if ((value3 < est_needed_step_err) and (found_h < 0.0)) then # if number 4
best_h = $glob_h
found_h = 1.0
end# end if 4
omniout_float(ALWAYS,"best_h",32,best_h,32,"")
opt_iter = opt_iter + 1
$glob_h = $glob_h * 0.5
end# end do number 2
if (found_h > 0.0) then # if number 4
$glob_h = best_h
else
omniout_str(ALWAYS,"No increment to obtain desired accuracy found")
end# end if 4
#END OPTIMIZE CODE
if ($glob_html_log) then # if number 4
html_log_file = File.new("html/entry.html","w")
end# end if 4
#BEGIN SOLUTION CODE
if (found_h > 0.0) then # if number 4
omniout_str(ALWAYS,"START of Soultion")
#Start Series -- INITIALIZE FOR SOLUTION
$array_x[1] = x_start
$array_x[2] = $glob_h
$glob_next_display = x_start
order_diff = 5
#Start Series $array_y2
term_no = 1
while (term_no <= order_diff) do # do number 2
$array_y2[term_no] = $array_y2_init[term_no] * expt($glob_h , (term_no - 1)) / factorial_1(term_no - 1)
term_no = term_no + 1
end# end do number 2
rows = order_diff
r_order = 1
while (r_order <= rows) do # do number 2
term_no = 1
while (term_no <= (rows - r_order + 1)) do # do number 3
it = term_no + r_order - 1
$array_y2_higher[r_order][term_no] = $array_y2_init[it]* expt($glob_h , (term_no - 1)) / ((factorial_1(term_no - 1)))
term_no = term_no + 1
end# end do number 3
r_order = r_order + 1
end# end do number 2
order_diff = 1
#Start Series $array_y1
term_no = 1
while (term_no <= order_diff) do # do number 2
$array_y1[term_no] = $array_y1_init[term_no] * expt($glob_h , (term_no - 1)) / factorial_1(term_no - 1)
term_no = term_no + 1
end# end do number 2
rows = order_diff
r_order = 1
while (r_order <= rows) do # do number 2
term_no = 1
while (term_no <= (rows - r_order + 1)) do # do number 3
it = term_no + r_order - 1
$array_y1_higher[r_order][term_no] = $array_y1_init[it]* expt($glob_h , (term_no - 1)) / ((factorial_1(term_no - 1)))
term_no = term_no + 1
end# end do number 3
r_order = r_order + 1
end# end do number 2
current_iter = 1
$glob_clock_start_sec = elapsed_time_seconds()
$glob_clock_sec = elapsed_time_seconds()
$glob_current_iter = 0
$glob_iter = 0
omniout_str(DEBUGL," ")
$glob_reached_optimal_h = true
$glob_optimal_clock_start_sec = elapsed_time_seconds()
while (($glob_current_iter < $glob_max_iter) and (($glob_check_sign * $array_x[1]) < ($glob_check_sign * x_end )) and ((convfloat($glob_clock_sec) - convfloat($glob_orig_start_sec)) < convfloat($glob_max_sec))) do # do number 2
#left paren 0001C
if (reached_interval()) then # if number 5
omniout_str(INFO," ")
omniout_str(INFO,"TOP MAIN SOLVE Loop")
end# end if 5
$glob_iter = $glob_iter + 1
$glob_clock_sec = elapsed_time_seconds()
$glob_current_iter = $glob_current_iter + 1
if ($glob_subiter_method == 1 ) then # if number 5
atomall()
elsif
($glob_subiter_method == 2 ) then # if number 6
subiter = 1
while (subiter <= 6) do # do number 3
atomall()
subiter = subiter + 1
end# end do number 3
else
subiter = 1
while (subiter <= 6 + $glob_max_terms) do # do number 3
atomall()
subiter = subiter + 1
end# end do number 3
end# end if 6
display_alot(current_iter)
if ($glob_look_poles) then # if number 6
#left paren 0004C
check_for_pole()
end# end if 6#was right paren 0004C
if (reached_interval()) then # if number 6
$glob_next_display = $glob_next_display + $glob_display_interval
end# end if 6
$array_x[1] = $array_x[1] + $glob_h
$array_x[2] = $glob_h
#Jump Series $array_y2
order_diff = 6
#START PART 1 SUM AND ADJUST
#START SUM AND ADJUST EQ =1
#sum_and_adjust $array_y2
#BEFORE ADJUST SUBSERIES EQ =1
ord = 6
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[6][iii] = $array_y2_higher[6][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 6
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 5
calc_term = 2
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[5][iii] = $array_y2_higher[5][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 5
calc_term = 2
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 5
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[5][iii] = $array_y2_higher[5][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 5
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 4
calc_term = 3
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[4][iii] = $array_y2_higher[4][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 4
calc_term = 3
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 4
calc_term = 2
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[4][iii] = $array_y2_higher[4][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 4
calc_term = 2
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 4
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[4][iii] = $array_y2_higher[4][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 4
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 3
calc_term = 4
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[3][iii] = $array_y2_higher[3][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 3
calc_term = 4
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 3
calc_term = 3
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[3][iii] = $array_y2_higher[3][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 3
calc_term = 3
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 3
calc_term = 2
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[3][iii] = $array_y2_higher[3][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 3
calc_term = 2
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 3
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[3][iii] = $array_y2_higher[3][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 3
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 5
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[2][iii] = $array_y2_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 2
calc_term = 5
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 4
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[2][iii] = $array_y2_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 2
calc_term = 4
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 3
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[2][iii] = $array_y2_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 2
calc_term = 3
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 2
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[2][iii] = $array_y2_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 2
calc_term = 2
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[2][iii] = $array_y2_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 2
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 6
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 6
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 5
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 5
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 4
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 4
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 3
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 3
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 2
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 2
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#BEFORE ADJUST SUBSERIES EQ =1
ord = 1
calc_term = 1
#adjust_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y2_higher_work[1][iii] = $array_y2_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =1
#BEFORE SUM SUBSERIES EQ =1
temp_sum = 0.0
ord = 1
calc_term = 1
#sum_subseries$array_y2
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y2_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y2_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =1
#END SUM AND ADJUST EQ =1
#END PART 1
#START PART 2 MOVE TERMS to REGULAR Array
term_no = $glob_max_terms
while (term_no >= 1) do # do number 3
$array_y2[term_no] = $array_y2_higher_work2[1][term_no]
ord = 1
while (ord <= order_diff) do # do number 4
$array_y2_higher[ord][term_no] = $array_y2_higher_work2[ord][term_no]
ord = ord + 1
end# end do number 4
term_no = term_no - 1
end# end do number 3
#END PART 2 HEVE MOVED TERMS to REGULAR Array
#Jump Series $array_y1
order_diff = 2
#START PART 1 SUM AND ADJUST
#START SUM AND ADJUST EQ =2
#sum_and_adjust $array_y1
#BEFORE ADJUST SUBSERIES EQ =2
ord = 2
calc_term = 1
#adjust_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y1_higher_work[2][iii] = $array_y1_higher[2][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =2
#BEFORE SUM SUBSERIES EQ =2
temp_sum = 0.0
ord = 2
calc_term = 1
#sum_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y1_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y1_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =2
#BEFORE ADJUST SUBSERIES EQ =2
ord = 1
calc_term = 2
#adjust_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y1_higher_work[1][iii] = $array_y1_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =2
#BEFORE SUM SUBSERIES EQ =2
temp_sum = 0.0
ord = 1
calc_term = 2
#sum_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y1_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y1_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =2
#BEFORE ADJUST SUBSERIES EQ =2
ord = 1
calc_term = 1
#adjust_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y1_higher_work[1][iii] = $array_y1_higher[1][iii] / expt($glob_h , (calc_term - 1)) / factorial_3(iii - calc_term , iii - 1)
iii = iii - 1
end# end do number 3
#AFTER ADJUST SUBSERIES EQ =2
#BEFORE SUM SUBSERIES EQ =2
temp_sum = 0.0
ord = 1
calc_term = 1
#sum_subseries$array_y1
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y1_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y1_higher_work2[ord][calc_term] = temp_sum * expt($glob_h , (calc_term - 1)) / (factorial_1(calc_term - 1))
#AFTER SUM SUBSERIES EQ =2
#END SUM AND ADJUST EQ =2
#END PART 1
#START PART 2 MOVE TERMS to REGULAR Array
term_no = $glob_max_terms
while (term_no >= 1) do # do number 3
$array_y1[term_no] = $array_y1_higher_work2[1][term_no]
ord = 1
while (ord <= order_diff) do # do number 4
$array_y1_higher[ord][term_no] = $array_y1_higher_work2[ord][term_no]
ord = ord + 1
end# end do number 4
term_no = term_no - 1
end# end do number 3
#END PART 2 HEVE MOVED TERMS to REGULAR Array
end# end do number 2#right paren 0001C
omniout_str(ALWAYS,"Finished!")
if ($glob_iter >= $glob_max_iter) then # if number 6
omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!")
end# end if 6
if (elapsed_time_seconds() - convfloat($glob_orig_start_sec) >= convfloat($glob_max_sec )) then # if number 6
omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!")
end# end if 6
$glob_clock_sec = elapsed_time_seconds()
omniout_str(INFO,"diff ( y2 , x , 5 ) = y1 ;")
omniout_str(INFO,"diff ( y1 , x , 1 ) = m1 * y2 ;")
omniout_int(INFO,"Iterations ",32,$glob_iter,4," ")
prog_report(x_start,x_end)
if ($glob_html_log) then # if number 6
logstart(html_log_file)
logitem_str(html_log_file,"2013-01-28T17:59:34-06:00")
logitem_str(html_log_file,"Ruby")
logitem_str(html_log_file,"mtest7")
logitem_str(html_log_file,"diff ( y2 , x , 5 ) = y1 ;")
logitem_float(html_log_file,x_start)
logitem_float(html_log_file,x_end)
logitem_float(html_log_file,$array_x[1])
logitem_float(html_log_file,$glob_h)
logitem_str(html_log_file,"16")
logitem_good_digits(html_log_file,$array_last_rel_error[1])
logitem_integer(html_log_file,$glob_max_terms)
logitem_float(html_log_file,$array_1st_rel_error[1])
logitem_float(html_log_file,$array_last_rel_error[1])
logitem_integer(html_log_file,$glob_iter)
logitem_pole(html_log_file,$array_type_pole[1])
if ($array_type_pole[1] == 1 or $array_type_pole[1] == 2) then # if number 7
logitem_float(html_log_file,$array_pole[1])
logitem_float(html_log_file,$array_pole[2])
0
else
logitem_str(html_log_file,"NA")
logitem_str(html_log_file,"NA")
0
end# end if 7
logitem_time(html_log_file,convfloat($glob_clock_sec))
if ($glob_percent_done < 100.0) then # if number 7
logitem_time(html_log_file,convfloat($glob_total_exp_sec))
0
else
logitem_str(html_log_file,"Done")
0
end# end if 7
log_revs(html_log_file," 165 | ")
logitem_str(html_log_file,"mtest7 diffeq.rb")
logitem_str(html_log_file,"mtest7 Ruby results")
logitem_str(html_log_file,"All Tests - All Languages")
logend(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logitem_str(html_log_file,"diff ( y1 , x , 1 ) = m1 * y2 ;")
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logitem_good_digits(html_log_file,$array_last_rel_error[2])
logditto(html_log_file)
logitem_float(html_log_file,$array_1st_rel_error[2])
logitem_float(html_log_file,$array_last_rel_error[2])
logditto(html_log_file)
logitem_pole(html_log_file,$array_type_pole[2])
if ($array_type_pole[2] == 1 or $array_type_pole[2] == 2) then # if number 7
logitem_float(html_log_file,$array_pole[1])
logitem_float(html_log_file,$array_pole[2])
0
else
logitem_str(html_log_file,"NA")
logitem_str(html_log_file,"NA")
0
end# end if 7
logditto(html_log_file)
if ($glob_percent_done < 100.0) then # if number 7
logditto(html_log_file)
0
else
logditto(html_log_file)
0
end# end if 7
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logditto(html_log_file)
logend(html_log_file)
end# end if 6
if ($glob_html_log) then # if number 6
html_log_file.close;
end# end if 6
end# end if 5
#END OUTFILEMAIN
end # End Function number 12
main()