#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) or (omniabs($array_y_higher[1][m-1]) < $glob_small_float) or (omniabs($array_y_higher[1][m-2]) < $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-1)*rm0-convfloat(m-2)*rm1 if (omniabs(hdrc) > $glob_small_float) then # if number 2 rcs = $glob_h/hdrc ord_no = convfloat(m-1)*rm0/hdrc - convfloat(m) + 2.0 $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 = false #TOP WHICH RADII EQ = 1 if ( not found 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 = true $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") end# end if 4 end# end if 3 end# end if 2 if ( not found 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] > 0.0) 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 = true $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") end# end if 4 end# end if 3 end# end if 2 if ( not found 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 = true $array_type_pole[1] = 3 if (reached_interval()) then # if number 3 omniout_str(ALWAYS,"NO POLE") end# end if 3 end# end if 2 if ( not found 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] > 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 = true $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") end# end if 4 end# end if 3 end# end if 2 if ( not found 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 = true if ($glob_display_flag) then # if number 3 if (reached_interval()) then # if number 4 omniout_str(ALWAYS,"Complex estimate of poles used") end# end if 4 end# end if 3 end# end if 2 if ( not found ) 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") 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 CONST - LINEAR $eq_no = 1 i = 1 $array_tmp1[1] = $array_const_0D1[1] * $array_x[1] #emit pre exp 1 $eq_no = 1 $array_tmp2[1] = exp($array_tmp1[1]) #emit pre mult CONST - LINEAR $eq_no = 1 i = 1 $array_tmp3[1] = $array_const_0D2[1] * $array_x[1] #emit pre exp 1 $eq_no = 1 $array_tmp4[1] = exp($array_tmp3[1]) #emit pre div FULL - FULL $eq_no = 1 i = 1 $array_tmp5[1] = ($array_tmp2[1] / ($array_tmp4[1])) #emit pre add CONST FULL $eq_no = 1 i = 1 $array_tmp6[1] = $array_const_0D0[1] + $array_tmp5[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_tmp6[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 CONST - LINEAR $eq_no = 1 i = 2 $array_tmp1[2] = $array_const_0D1[1] * $array_x[2] #emit pre exp ID_LINEAR iii = 2 $eq_no = 1 $array_tmp2[2] = $array_tmp2[1] * $array_tmp1[2] / 1 #emit pre mult CONST - LINEAR $eq_no = 1 i = 2 $array_tmp3[2] = $array_const_0D2[1] * $array_x[2] #emit pre exp ID_LINEAR iii = 2 $eq_no = 1 $array_tmp4[2] = $array_tmp4[1] * $array_tmp3[2] / 1 #emit pre div FULL - FULL $eq_no = 1 i = 2 $array_tmp5[2] = (($array_tmp2[2] - ats(2,$array_tmp4,$array_tmp5,2))/$array_tmp4[1]) #emit pre add CONST FULL $eq_no = 1 i = 2 $array_tmp6[2] = $array_tmp5[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_tmp6[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 exp ID_LINEAR iii = 3 $eq_no = 1 $array_tmp2[3] = $array_tmp2[2] * $array_tmp1[2] / 2 #emit pre exp ID_LINEAR iii = 3 $eq_no = 1 $array_tmp4[3] = $array_tmp4[2] * $array_tmp3[2] / 2 #emit pre div FULL - FULL $eq_no = 1 i = 3 $array_tmp5[3] = (($array_tmp2[3] - ats(3,$array_tmp4,$array_tmp5,2))/$array_tmp4[1]) #emit pre add CONST FULL $eq_no = 1 i = 3 $array_tmp6[3] = $array_tmp5[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_tmp6[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 exp ID_LINEAR iii = 4 $eq_no = 1 $array_tmp2[4] = $array_tmp2[3] * $array_tmp1[2] / 3 #emit pre exp ID_LINEAR iii = 4 $eq_no = 1 $array_tmp4[4] = $array_tmp4[3] * $array_tmp3[2] / 3 #emit pre div FULL - FULL $eq_no = 1 i = 4 $array_tmp5[4] = (($array_tmp2[4] - ats(4,$array_tmp4,$array_tmp5,2))/$array_tmp4[1]) #emit pre add CONST FULL $eq_no = 1 i = 4 $array_tmp6[4] = $array_tmp5[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_tmp6[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 exp ID_LINEAR iii = 5 $eq_no = 1 $array_tmp2[5] = $array_tmp2[4] * $array_tmp1[2] / 4 #emit pre exp ID_LINEAR iii = 5 $eq_no = 1 $array_tmp4[5] = $array_tmp4[4] * $array_tmp3[2] / 4 #emit pre div FULL - FULL $eq_no = 1 i = 5 $array_tmp5[5] = (($array_tmp2[5] - ats(5,$array_tmp4,$array_tmp5,2))/$array_tmp4[1]) #emit pre add CONST FULL $eq_no = 1 i = 5 $array_tmp6[5] = $array_tmp5[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_tmp6[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 exp LINEAR $eq_no = 1 $array_tmp2[kkk] = $array_tmp2[kkk - 1] * $array_tmp1[2] / (kkk - 1) #emit exp LINEAR $eq_no = 1 $array_tmp4[kkk] = $array_tmp4[kkk - 1] * $array_tmp3[2] / (kkk - 1) #emit div FULL FULL $eq_no = 1 $array_tmp5[kkk] = (($array_tmp2[kkk] - ats(kkk,$array_tmp4,$array_tmp5,2))/$array_tmp4[1]) #emit NOT FULL - FULL add $eq_no = 1 $array_tmp6[kkk] = $array_tmp5[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_tmp6[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() 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 15 # Begin Function number 16 def logditto(file) file.printf("") file.printf("ditto") file.printf("") end # End Function number 16 # Begin Function number 17 def logitem_integer(file,n) file.printf("") file.printf("%d",n) file.printf("") end # End Function number 17 # Begin Function number 18 def logitem_str(file,str) file.printf("") file.printf(str) file.printf("") end # End Function number 18 # Begin Function number 19 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 19 # Begin Function number 20 def log_revs(file,revs) file.printf(revs.to_s) end # End Function number 20 # Begin Function number 21 def logitem_float(file,x) file.printf("") file.printf(x.to_s) file.printf("") end # End Function number 21 # Begin Function number 22 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 22 # Begin Function number 23 def logstart(file) file.printf("") end # End Function number 23 # Begin Function number 24 def logend(file) file.printf("") end # End Function number 24 # Begin Function number 25 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 25 # Begin Function number 26 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 26 # Begin Function number 27 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 27 # Begin Function number 28 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 28 # Begin Function number 29 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 29 # Begin Function number 30 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 30 # Begin Function number 31 def convfp(mmm) return(mmm) end # End Function number 31 # Begin Function number 32 def convfloat(mmm) return(mmm) end # End Function number 32 # Begin Function number 33 def elapsed_time_seconds() t = Time.now return(t.to_i) end # End Function number 33 # Begin Function number 34 def expt(x,y) return(x**y) end # End Function number 34 # Begin Function number 35 def Si(x) return(0.0) end # End Function number 35 # Begin Function number 36 def Ci(x) return(0.0) end # End Function number 36 # Begin Function number 37 def abs(x) return(x.abs) end # End Function number 37 # Begin Function number 38 def omniabs(x) return(x.abs) end # End Function number 38 # Begin Function number 39 def exp(x) return(Math.exp(x)) end # End Function number 39 # Begin Function number 40 def trunc(x) return(x.to_i) end # End Function number 40 # Begin Function number 41 def floor(x) return(x.floor) end # End Function number 41 # Begin Function number 42 def sin(x) return(Math.sin(x)) end # End Function number 42 # Begin Function number 43 def cos(x) return(Math.cos(x)) end # End Function number 43 # Begin Function number 44 def tan(x) return(Math.tan(x)) end # End Function number 44 # Begin Function number 45 def arccos(x) return(Math.acos(x)) end # End Function number 45 # Begin Function number 46 def arccosh(x) return(Math.acosh(x)) end # End Function number 46 # Begin Function number 47 def arcsin(x) return(Math.asin(x)) end # End Function number 47 # Begin Function number 48 def arcsinh(x) return(Math.asinh(x)) end # End Function number 48 # Begin Function number 49 def arctan(x) return(Math.atan(x)) end # End Function number 49 # Begin Function number 50 def arctanh(x) return(Math.atanh(x)) end # End Function number 50 # Begin Function number 51 def cosh(x) return(Math.cosh(x)) end # End Function number 51 # Begin Function number 52 def erf(x) return(Math.erf(x)) end # End Function number 52 # Begin Function number 53 def ln(x) return(Math.log(x)) end # End Function number 53 # Begin Function number 54 def log10(x) return(Math.log10(x)) end # End Function number 54 # Begin Function number 55 def sinh(x) return(Math.sinh(x)) end # End Function number 55 # Begin Function number 56 def tanh(x) return(Math.tanh(x)) end # End Function number 56 # Begin Function number 57 def sqrt(x) return(Math.sqrt(x)) end # End Function number 57 # Begin Function number 58 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 58 # Begin Function number 59 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 59 #END ATS LIBRARY BLOCK #BEGIN USER DEF BLOCK #BEGIN USER DEF BLOCK def exact_soln_y (x) return(- 10.0 * (exp(0.1* x)/exp(0.2*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_log10normmin = 0.1 $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_hmax = 1.0 $glob_hmin = 0.00000000001 $glob_hmin_init = 0.001 $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_log10_abserr = 0.1e-10 $glob_log10_relerr = 0.1e-10 $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-50 $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_log10abserr = 0.0 $glob_log10relerr = 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/div_exp_exppostode.ode#################") omniout_str(ALWAYS,"diff ( y , x , 1 ) = exp(0.1 * x) / exp(0.2 * x);") 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=-5.0 ") omniout_str(ALWAYS,"x_end=5.0 ") omniout_str(ALWAYS,"$array_y_init[0 + 1] = exact_soln_y(x_start)") omniout_str(ALWAYS,"$glob_h=0.05 ") omniout_str(ALWAYS,"$glob_look_poles=true ") omniout_str(ALWAYS,"$glob_max_iter=1000000 ") omniout_str(ALWAYS,"$glob_display_interval=0.1 ") omniout_str(ALWAYS,"$glob_max_minutes=10 ") 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(- 10.0 * (exp(0.1* x)/exp(0.2*x))) ") omniout_str(ALWAYS,"end") 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 $glob_log10_abserr = -8.0 $glob_log10_relerr = -8.0 $glob_hmax = 0.01 #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_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_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_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_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_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_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_0D1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 2 $array_const_0D1[term] = 0.0 term = term + 1 end# end do number 2 $array_const_0D1[1] = 0.1 $array_const_0D2 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 2 $array_const_0D2[term] = 0.0 term = term + 1 end# end do number 2 $array_const_0D2[1] = 0.2 $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=-5.0 x_end=5.0 $array_y_init[0 + 1] = exact_soln_y(x_start) $glob_h=0.05 $glob_look_poles=true $glob_max_iter=1000000 $glob_display_interval=0.1 $glob_max_minutes=10 #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) $glob_abserr = expt(10.0 , ($glob_log10_abserr)) $glob_relerr = expt(10.0 , ($glob_log10_relerr)) 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 #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 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_log10normmin = -$glob_large_float if (omniabs($array_y_higher[1][1]) > $glob_small_float) then # if number 3 tmp = omniabs($array_y_higher[1][1]) log10norm = (log10(tmp)) if (log10norm < $glob_log10normmin) then # if number 4 $glob_log10normmin = log10norm end# end if 4 end# end if 3 display_alot(current_iter) $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 ) = exp(0.1 * x) / exp(0.2 * x);") 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-12T21:52:22-06:00") logitem_str(html_log_file,"Ruby") logitem_str(html_log_file,"div_exp_exp") logitem_str(html_log_file,"diff ( y , x , 1 ) = exp(0.1 * x) / exp(0.2 * x);") 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," 156 ") logitem_str(html_log_file,"div_exp_exp diffeq.rb") logitem_str(html_log_file,"div_exp_exp Ruby results") logitem_str(html_log_file,"Languages compared - single equations") 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()