#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_y[1]) < min_size) then # if number 1
min_size = omniabs($array_y[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_y[no_terms-3] + $array_y[no_terms - 2] * hn_div_ho + $array_y[no_terms - 1] * hn_div_ho_2 + $array_y[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_y(ind_var)
omniout_float(ALWAYS,"y[1] (analytic) ",33,analytic_val_y,20," ")
term_no = 1
numeric_val = $array_y[term_no]
abserr = omniabs(numeric_val - analytic_val_y)
omniout_float(ALWAYS,"y[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," ")
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_y_higher[1][1]) > $glob_small_float) then # if number 1
tmp = omniabs($array_y_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 - 1 - 1
while ((m >= 10) and ((omniabs($array_y_higher[1][m]) < $glob_small_float * $glob_small_float) or (omniabs($array_y_higher[1][m-1]) < $glob_small_float * $glob_small_float) or (omniabs($array_y_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_y_higher[1][m]/$array_y_higher[1][m-1]
rm1 = $array_y_higher[1][m-1]/$array_y_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
#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 - 1 - 1
cnt = 0
while ((cnt < 5) and (n >= 10)) do # do number 2
if (omniabs($array_y_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_y_higher[1][m]) >= ($glob_large_float)) or (omniabs($array_y_higher[1][m-1]) >=($glob_large_float)) or (omniabs($array_y_higher[1][m-2]) >= ($glob_large_float)) or (omniabs($array_y_higher[1][m-3]) >= ($glob_large_float)) or (omniabs($array_y_higher[1][m-4]) >= ($glob_large_float)) or (omniabs($array_y_higher[1][m-5]) >= ($glob_large_float))) or ((omniabs($array_y_higher[1][m]) <= ($glob_small_float)) or (omniabs($array_y_higher[1][m-1]) <=($glob_small_float)) or (omniabs($array_y_higher[1][m-2]) <= ($glob_small_float)) or (omniabs($array_y_higher[1][m-3]) <= ($glob_small_float)) or (omniabs($array_y_higher[1][m-4]) <= ($glob_small_float)) or (omniabs($array_y_higher[1][m-5]) <= ($glob_small_float)))) then # if number 2
rad_c = $glob_large_float
ord_no = $glob_large_float
else
rm0 = ($array_y_higher[1][m])/($array_y_higher[1][m-1])
rm1 = ($array_y_higher[1][m-1])/($array_y_higher[1][m-2])
rm2 = ($array_y_higher[1][m-2])/($array_y_higher[1][m-3])
rm3 = ($array_y_higher[1][m-3])/($array_y_higher[1][m-4])
rm4 = ($array_y_higher[1][m-4])/($array_y_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
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 2
$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 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"Complex estimate of poles used for equation 1")
end# end if 4
end# end if 3
end# end if 2
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 2
$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 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"Real estimate of pole used for equation 1")
end# end if 4
end# end if 3
end# end if 2
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 2
$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 3
omniout_str(ALWAYS,"NO POLE for equation 1")
end# end if 3
end# end if 2
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 2
$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 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"Real estimate of pole used for equation 1")
end# end if 4
end# end if 3
end# end if 2
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 2
$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 3
if (reached_interval()) then # if number 4
omniout_str(ALWAYS,"Complex estimate of poles used for equation 1")
end# end if 4
end# end if 3
end# end if 2
if (1 != found_sing ) then # if number 2
$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 3
omniout_str(ALWAYS,"NO POLE for equation 1")
end# end if 3
end# end if 2
#BOTTOM WHICH RADII EQ = 1
$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 2
$array_pole[1] = $array_poles[1][1]
$array_pole[2] = $array_poles[1][2]
end# end if 2
#BOTTOM WHICH RADIUS EQ = 1
#START ADJUST ALL SERIES
if ($array_pole[1] * $glob_ratio_of_radius < omniabs($glob_h)) then # if number 2
h_new = $array_pole[1] * $glob_ratio_of_radius
term = 1
ratio = 1.0
while (term <= $glob_max_terms) do # do number 2
$array_y[term] = $array_y[term]* ratio
$array_y_higher[1][term] = $array_y_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 2
#BOTTOM ADJUST ALL SERIES
if (reached_interval()) then # if number 2
display_pole()
end# end if 2
end # End Function number 10
# Begin Function number 11
def get_norms()
if ( not $glob_initial_pass) then # if number 2
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_y[iii]) > $array_norms[iii]) then # if number 3
$array_norms[iii] = omniabs($array_y[iii])
end# end if 3
iii = iii + 1
end# end do number 2
#BOTTOM GET NORMS
end# end if 2
end # End Function number 11
# Begin Function number 12
def atomall()
#TOP ATOMALL
#END OUTFILE1
#BEGIN ATOMHDR1
#emit pre mult FULL CONST $eq_no = 1 i = 1
$array_tmp1[1] = $array_m1[1] * $array_const_2D0[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 1
$array_tmp2[1] = $array_x[1] - $array_const_6D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 1
$array_tmp3[1] = $array_tmp1[1] / $array_tmp2[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 1
$array_tmp4[1] = $array_x[1] - $array_const_6D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 1
$array_tmp5[1] = $array_tmp3[1] / $array_tmp4[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 1
$array_tmp6[1] = $array_x[1] - $array_const_6D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 1
$array_tmp7[1] = $array_tmp5[1] / $array_tmp6[1]
#emit pre add CONST FULL $eq_no = 1 i = 1
$array_tmp8[1] = $array_const_0D0[1] + $array_tmp7[1]
#emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5
if ( not $array_y_set_initial[1][2]) then # if number 1
if (1 <= $glob_max_terms) then # if number 2
temporary = $array_tmp8[1] * expt($glob_h , (1)) * factorial_3(0,1)
$array_y[2] = temporary
$array_y_higher[1][2] = temporary
temporary = temporary / $glob_h * (1.0)
$array_y_higher[2][1] = temporary
end# end if 2
end# end if 1
kkk = 2
#END ATOMHDR1
#BEGIN ATOMHDR2
#emit pre mult FULL CONST $eq_no = 1 i = 2
$array_tmp1[2] = $array_m1[2] * $array_const_2D0[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 2
$array_tmp2[2] = $array_x[2]
#emit pre div FULL - LINEAR $eq_no = 1 i = 2
$array_tmp3[2] = ($array_tmp1[2] - $array_tmp3[1] * $array_tmp2[2]) / $array_tmp2[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 2
$array_tmp4[2] = $array_x[2]
#emit pre div FULL - LINEAR $eq_no = 1 i = 2
$array_tmp5[2] = ($array_tmp3[2] - $array_tmp5[1] * $array_tmp4[2]) / $array_tmp4[1]
#emit pre sub LINEAR - CONST $eq_no = 1 i = 2
$array_tmp6[2] = $array_x[2]
#emit pre div FULL - LINEAR $eq_no = 1 i = 2
$array_tmp7[2] = ($array_tmp5[2] - $array_tmp7[1] * $array_tmp6[2]) / $array_tmp6[1]
#emit pre add CONST FULL $eq_no = 1 i = 2
$array_tmp8[2] = $array_tmp7[2]
#emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5
if ( not $array_y_set_initial[1][3]) then # if number 1
if (2 <= $glob_max_terms) then # if number 2
temporary = $array_tmp8[2] * expt($glob_h , (1)) * factorial_3(1,2)
$array_y[3] = temporary
$array_y_higher[1][3] = temporary
temporary = temporary / $glob_h * (2.0)
$array_y_higher[2][2] = temporary
end# end if 2
end# end if 1
kkk = 3
#END ATOMHDR2
#BEGIN ATOMHDR3
#emit pre mult FULL CONST $eq_no = 1 i = 3
$array_tmp1[3] = $array_m1[3] * $array_const_2D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 3
$array_tmp3[3] = ($array_tmp1[3] - $array_tmp3[2] * $array_tmp2[2]) / $array_tmp2[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 3
$array_tmp5[3] = ($array_tmp3[3] - $array_tmp5[2] * $array_tmp4[2]) / $array_tmp4[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 3
$array_tmp7[3] = ($array_tmp5[3] - $array_tmp7[2] * $array_tmp6[2]) / $array_tmp6[1]
#emit pre add CONST FULL $eq_no = 1 i = 3
$array_tmp8[3] = $array_tmp7[3]
#emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5
if ( not $array_y_set_initial[1][4]) then # if number 1
if (3 <= $glob_max_terms) then # if number 2
temporary = $array_tmp8[3] * expt($glob_h , (1)) * factorial_3(2,3)
$array_y[4] = temporary
$array_y_higher[1][4] = temporary
temporary = temporary / $glob_h * (3.0)
$array_y_higher[2][3] = temporary
end# end if 2
end# end if 1
kkk = 4
#END ATOMHDR3
#BEGIN ATOMHDR4
#emit pre mult FULL CONST $eq_no = 1 i = 4
$array_tmp1[4] = $array_m1[4] * $array_const_2D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 4
$array_tmp3[4] = ($array_tmp1[4] - $array_tmp3[3] * $array_tmp2[2]) / $array_tmp2[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 4
$array_tmp5[4] = ($array_tmp3[4] - $array_tmp5[3] * $array_tmp4[2]) / $array_tmp4[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 4
$array_tmp7[4] = ($array_tmp5[4] - $array_tmp7[3] * $array_tmp6[2]) / $array_tmp6[1]
#emit pre add CONST FULL $eq_no = 1 i = 4
$array_tmp8[4] = $array_tmp7[4]
#emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5
if ( not $array_y_set_initial[1][5]) then # if number 1
if (4 <= $glob_max_terms) then # if number 2
temporary = $array_tmp8[4] * expt($glob_h , (1)) * factorial_3(3,4)
$array_y[5] = temporary
$array_y_higher[1][5] = temporary
temporary = temporary / $glob_h * (4.0)
$array_y_higher[2][4] = temporary
end# end if 2
end# end if 1
kkk = 5
#END ATOMHDR4
#BEGIN ATOMHDR5
#emit pre mult FULL CONST $eq_no = 1 i = 5
$array_tmp1[5] = $array_m1[5] * $array_const_2D0[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 5
$array_tmp3[5] = ($array_tmp1[5] - $array_tmp3[4] * $array_tmp2[2]) / $array_tmp2[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 5
$array_tmp5[5] = ($array_tmp3[5] - $array_tmp5[4] * $array_tmp4[2]) / $array_tmp4[1]
#emit pre div FULL - LINEAR $eq_no = 1 i = 5
$array_tmp7[5] = ($array_tmp5[5] - $array_tmp7[4] * $array_tmp6[2]) / $array_tmp6[1]
#emit pre add CONST FULL $eq_no = 1 i = 5
$array_tmp8[5] = $array_tmp7[5]
#emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5
if ( not $array_y_set_initial[1][6]) then # if number 1
if (5 <= $glob_max_terms) then # if number 2
temporary = $array_tmp8[5] * expt($glob_h , (1)) * factorial_3(4,5)
$array_y[6] = temporary
$array_y_higher[1][6] = temporary
temporary = temporary / $glob_h * (5.0)
$array_y_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 mult FULL CONST $eq_no = 1 i = 1
$array_tmp1[kkk] = $array_m1[kkk] * $array_const_2D0[1]
#emit div FULL LINEAR $eq_no = 1 i = 1
$array_tmp3[kkk] = -ats(kkk,$array_tmp2,$array_tmp3,2) / $array_tmp2[1]
#emit div FULL LINEAR $eq_no = 1 i = 1
$array_tmp5[kkk] = -ats(kkk,$array_tmp4,$array_tmp5,2) / $array_tmp4[1]
#emit div FULL LINEAR $eq_no = 1 i = 1
$array_tmp7[kkk] = -ats(kkk,$array_tmp6,$array_tmp7,2) / $array_tmp6[1]
#emit NOT FULL - FULL add $eq_no = 1
$array_tmp8[kkk] = $array_tmp7[kkk]
#emit assign $eq_no = 1
order_d = 1
if (kkk + order_d + 1 <= $glob_max_terms) then # if number 1
if ( not $array_y_set_initial[1][kkk + order_d]) then # if number 2
temporary = $array_tmp8[kkk] * expt($glob_h , (order_d)) * factorial_3((kkk - 1),(kkk + order_d - 1))
$array_y[kkk + order_d] = temporary
$array_y_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_y_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_y (x)
return(1.0/ ( x - 6.0 ) / ( x - 6.0 ))
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 = 1
$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/sing6postode.ode#################")
omniout_str(ALWAYS,"diff ( y , x , 1 ) = m1 * 2.0 / ( x - 6.0 ) / ( x - 6.0 ) / ( x - 6.0) ;")
omniout_str(ALWAYS,"!")
omniout_str(ALWAYS,"#BEGIN FIRST INPUT BLOCK")
omniout_str(ALWAYS,"# Digits:=64; ELIMINATED in preodein.rb")
omniout_str(ALWAYS,"max_terms=40 ")
omniout_str(ALWAYS,"!")
omniout_str(ALWAYS,"#END FIRST INPUT BLOCK")
omniout_str(ALWAYS,"#BEGIN SECOND INPUT BLOCK")
omniout_str(ALWAYS,"x_start=2.0 ")
omniout_str(ALWAYS,"x_end=3.0 ")
omniout_str(ALWAYS,"$array_y_init[0 + 1] = exact_soln_y(x_start)")
omniout_str(ALWAYS,"$glob_look_poles=true ")
omniout_str(ALWAYS,"$glob_max_iter=100 ")
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_y (x)")
omniout_str(ALWAYS,"return(1.0/ ( x - 6.0 ) / ( x - 6.0 )) ")
omniout_str(ALWAYS,"end")
omniout_str(ALWAYS,"")
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:=64; ELIMINATED in preodein.rb
max_terms=40
#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_y_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_y= Array.new(max_terms + 1)
$array_x= 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_tmp4= Array.new(max_terms + 1)
$array_tmp5= Array.new(max_terms + 1)
$array_tmp6= Array.new(max_terms + 1)
$array_tmp7= Array.new(max_terms + 1)
$array_tmp8= Array.new(max_terms + 1)
$array_m1= Array.new(max_terms + 1)
$array_y_higher = array2d(2 + 1 ,max_terms + 1 )
$array_y_higher_work = array2d(2 + 1 ,max_terms + 1 )
$array_y_higher_work2 = array2d(2 + 1 ,max_terms + 1 )
$array_y_set_initial = array2d(2 + 1 ,max_terms + 1 )
$array_poles = array2d(1 + 1 ,3 + 1 )
$array_real_pole = array2d(1 + 1 ,3 + 1 )
$array_complex_pole = array2d(1 + 1 ,3 + 1 )
$array_fact_2 = array2d(max_terms + 1 ,max_terms + 1 )
term = 1
while (term <= max_terms) do # do number 2
$array_y_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_y[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_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_tmp4[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp5[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp6[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp7[term] = 0.0
term = term + 1
end# end do number 2
term = 1
while (term <= max_terms) do # do number 2
$array_tmp8[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 <=2) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y_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_y_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_y_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 <=2) do # do number 2
term = 1
while (term <= max_terms) do # do number 3
$array_y_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 <=1) 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 <=1) 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 <=1) 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_y = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_y[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_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_tmp4 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp4[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp5 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp5[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp6 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp6[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp7 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp7[term] = 0.0
term = term + 1
end# end do number 2
$array_tmp8 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_tmp8[term] = 0.0
term = term + 1
end# end do number 2
$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_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_2D0 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_const_2D0[term] = 0.0
term = term + 1
end# end do number 2
$array_const_2D0[1] = 2.0
$array_const_6D0 = Array.new(max_terms+1 + 1)
term = 1
while (term <= max_terms + 1) do # do number 2
$array_const_6D0[term] = 0.0
term = term + 1
end# end do number 2
$array_const_6D0[1] = 6.0
$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=2.0
x_end=3.0
$array_y_init[0 + 1] = exact_soln_y(x_start)
$glob_look_poles=true
$glob_max_iter=100
#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_y_set_initial[1][1] = true
$array_y_set_initial[1][2] = false
$array_y_set_initial[1][3] = false
$array_y_set_initial[1][4] = false
$array_y_set_initial[1][5] = false
$array_y_set_initial[1][6] = false
$array_y_set_initial[1][7] = false
$array_y_set_initial[1][8] = false
$array_y_set_initial[1][9] = false
$array_y_set_initial[1][10] = false
$array_y_set_initial[1][11] = false
$array_y_set_initial[1][12] = false
$array_y_set_initial[1][13] = false
$array_y_set_initial[1][14] = false
$array_y_set_initial[1][15] = false
$array_y_set_initial[1][16] = false
$array_y_set_initial[1][17] = false
$array_y_set_initial[1][18] = false
$array_y_set_initial[1][19] = false
$array_y_set_initial[1][20] = false
$array_y_set_initial[1][21] = false
$array_y_set_initial[1][22] = false
$array_y_set_initial[1][23] = false
$array_y_set_initial[1][24] = false
$array_y_set_initial[1][25] = false
$array_y_set_initial[1][26] = false
$array_y_set_initial[1][27] = false
$array_y_set_initial[1][28] = false
$array_y_set_initial[1][29] = false
$array_y_set_initial[1][30] = false
$array_y_set_initial[1][31] = false
$array_y_set_initial[1][32] = false
$array_y_set_initial[1][33] = false
$array_y_set_initial[1][34] = false
$array_y_set_initial[1][35] = false
$array_y_set_initial[1][36] = false
$array_y_set_initial[1][37] = false
$array_y_set_initial[1][38] = false
$array_y_set_initial[1][39] = false
$array_y_set_initial[1][40] = 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 2
$glob_h = $glob_display_interval
end# end if 2
if ($glob_max_h < $glob_h) then # if number 2
$glob_h = $glob_max_h
end# end if 2
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 = 1
#Start Series $array_y
term_no = 1
while (term_no <= order_diff) do # do number 3
$array_y[term_no] = $array_y_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_y_higher[r_order][term_no] = $array_y_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
atomall()
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 2
best_h = $glob_h
found_h = 1.0
end# end if 2
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 2
$glob_h = best_h
else
omniout_str(ALWAYS,"No increment to obtain desired accuracy found")
end# end if 2
#END OPTIMIZE CODE
if ($glob_html_log) then # if number 2
html_log_file = File.new("html/entry.html","w")
end# end if 2
#BEGIN SOLUTION CODE
if (found_h > 0.0) then # if number 2
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 = 1
#Start Series $array_y
term_no = 1
while (term_no <= order_diff) do # do number 2
$array_y[term_no] = $array_y_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_y_higher[r_order][term_no] = $array_y_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 3
omniout_str(INFO," ")
omniout_str(INFO,"TOP MAIN SOLVE Loop")
end# end if 3
$glob_iter = $glob_iter + 1
$glob_clock_sec = elapsed_time_seconds()
$glob_current_iter = $glob_current_iter + 1
atomall()
display_alot(current_iter)
if ($glob_look_poles) then # if number 3
#left paren 0004C
check_for_pole()
end# end if 3#was right paren 0004C
if (reached_interval()) then # if number 3
$glob_next_display = $glob_next_display + $glob_display_interval
end# end if 3
$array_x[1] = $array_x[1] + $glob_h
$array_x[2] = $glob_h
#Jump Series $array_y
order_diff = 2
#START PART 1 SUM AND ADJUST
#START SUM AND ADJUST EQ =1
#sum_and_adjust $array_y
#BEFORE ADJUST SUBSERIES EQ =1
ord = 2
calc_term = 1
#adjust_subseries$array_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y_higher_work[2][iii] = $array_y_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_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y_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_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y_higher_work[1][iii] = $array_y_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_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y_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_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
$array_y_higher_work[1][iii] = $array_y_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_y
iii = $glob_max_terms
while (iii >= calc_term) do # do number 3
temp_sum = temp_sum + $array_y_higher_work[ord][iii]
iii = iii - 1
end# end do number 3
$array_y_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_y[term_no] = $array_y_higher_work2[1][term_no]
ord = 1
while (ord <= order_diff) do # do number 4
$array_y_higher[ord][term_no] = $array_y_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 3
omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!")
end# end if 3
if (elapsed_time_seconds() - convfloat($glob_orig_start_sec) >= convfloat($glob_max_sec )) then # if number 3
omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!")
end# end if 3
$glob_clock_sec = elapsed_time_seconds()
omniout_str(INFO,"diff ( y , x , 1 ) = m1 * 2.0 / ( x - 6.0 ) / ( x - 6.0 ) / ( x - 6.0) ;")
omniout_int(INFO,"Iterations ",32,$glob_iter,4," ")
prog_report(x_start,x_end)
if ($glob_html_log) then # if number 3
logstart(html_log_file)
logitem_str(html_log_file,"2013-01-28T19:21:40-06:00")
logitem_str(html_log_file,"Ruby")
logitem_str(html_log_file,"sing6")
logitem_str(html_log_file,"diff ( y , x , 1 ) = m1 * 2.0 / ( x - 6.0 ) / ( x - 6.0 ) / ( x - 6.0) ;")
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 4
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 4
logitem_time(html_log_file,convfloat($glob_clock_sec))
if ($glob_percent_done < 100.0) then # if number 4
logitem_time(html_log_file,convfloat($glob_total_exp_sec))
0
else
logitem_str(html_log_file,"Done")
0
end# end if 4
log_revs(html_log_file," 165 | ")
logitem_str(html_log_file,"sing6 diffeq.rb")
logitem_str(html_log_file,"sing6 Ruby results")
logitem_str(html_log_file,"All Tests - All Languages")
logend(html_log_file)
end# end if 3
if ($glob_html_log) then # if number 3
html_log_file.close;
end# end if 3
end# end if 2
#END OUTFILEMAIN
end # End Function number 12
main()