#BEGIN OUTFILE1 ALWAYS = 1 INFO = 2 DEBUGL = 3 DEBUGMASSIVE = 4 MAX_UNCHANGED = 10 # Begin Function number 3 def display_poles() if ($glob_type_given_pole == 4) then # if number 1 rad_given = sqrt(expt($array_x[1] - $array_given_rad_poles[1][1],2.0) + expt($array_given_rad_poles[1][2],2.0)) omniout_float(ALWAYS,"Radius of convergence (given) for eq 1 ",4, rad_given,4," ") omniout_float(ALWAYS,"Order of pole (given) ",4, $array_given_ord_poles[1][1],4," ") elsif ($glob_type_given_pole == 3) then # if number 2 omniout_str(ALWAYS,"NO POLE (given) for Equation 1") else omniout_str(ALWAYS,"NO INFO (given) for Equation 1") end# end if 2 if ($array_poles[1][1] != $glob_large_float) then # if number 2 omniout_float(ALWAYS,"Radius of convergence (ratio test) for eq 1 ",4, $array_poles[1][1],4," ") omniout_str(ALWAYS,"Order of pole (ratio test) Not computed") else omniout_str(ALWAYS,"NO POLE (ratio test) for Equation 1") end# end if 2 if (($array_real_poles[1][1] > 0.0) and ($array_real_poles[1][1] != $glob_large_float)) then # if number 2 omniout_float(ALWAYS,"Radius of convergence (three term test) for eq 1 ",4, $array_real_poles[1][1],4," ") omniout_float(ALWAYS,"Order of pole (three term test) ",4, $array_real_poles[1][2],4," ") else omniout_str(ALWAYS,"NO REAL POLE (three term test) for Equation 1") end# end if 2 if (($array_complex_poles[1][1] > 0.0) and ($array_complex_poles[1][1] != $glob_large_float)) then # if number 2 omniout_float(ALWAYS,"Radius of convergence (six term test) for eq 1 ",4, $array_complex_poles[1][1],4," ") omniout_float(ALWAYS,"Order of pole (six term test) ",4, $array_complex_poles[1][2],4," ") else omniout_str(ALWAYS,"NO COMPLEX POLE (six term test) for Equation 1") end# end if 2 end # End Function number 3 # Begin Function number 4 def check_sign( x0 ,xf) if (xf > x0) then # if number 2 ret = 1.0 else ret = -1.0 end# end if 2 return ( ret) end # End Function number 4 # Begin Function number 5 def est_size_answer() min_size = $glob_large_float if (omniabs($array_y[1]) < min_size) then # if number 2 min_size = omniabs($array_y[1]) omniout_float(ALWAYS,"min_size",32,min_size,32,"") end# end if 2 if (min_size < 1.0) then # if number 2 min_size = 1.0 omniout_float(ALWAYS,"min_size",32,min_size,32,"") end# end if 2 return ( min_size) end # End Function number 5 # Begin Function number 6 def test_suggested_h() max_estimated_step_error = 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,"") est_tmp = 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 (est_tmp >= max_estimated_step_error) then # if number 2 max_estimated_step_error = est_tmp end# end if 2 omniout_float(ALWAYS,"max_estimated_step_error",32,max_estimated_step_error,32,"") return ( max_estimated_step_error) end # End Function number 6 # Begin Function number 7 def reached_interval() if ($glob_check_sign * ($array_x[1]) >= $glob_check_sign * $glob_next_display) then # if number 2 ret = true else ret = false end# end if 2 return(ret) end # End Function number 7 # Begin Function number 8 def display_alot(iter) #TOP DISPLAY ALOT if (reached_interval()) then # if number 2 if (iter >= 0) then # if number 3 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 4 relerr = abserr*100.0/omniabs(analytic_val_y) if (relerr > 0.0000000000000000000000000000000001) then # if number 5 $glob_good_digits = -trunc(log10(relerr)) + 3 else $glob_good_digits = 16 end# end if 5 else relerr = -1.0 $glob_good_digits = -1 end# end if 4 if ($glob_iter == 1) then # if number 4 $array_1st_rel_error[1] = relerr else $array_last_rel_error[1] = relerr end# end if 4 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 3 #BOTTOM DISPLAY ALOT end# end if 2 end # End Function number 8 # Begin Function number 9 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 2 tmp = omniabs($array_y_higher[1][1]) if (tmp < $glob_normmax) then # if number 3 $glob_normmax = tmp end# end if 3 end# end if 2 if ($glob_look_poles and (omniabs($array_pole[1]) > $glob_small_float) and ($array_pole[1] != $glob_large_float)) then # if number 2 sz2 = $array_pole[1]/10.0 if (sz2 < hnew) then # if number 3 omniout_float(INFO,"$glob_h adjusted to ",20,h_param,12,"due to singularity.") omniout_str(INFO,"Reached Optimal") return(hnew) end# end if 3 end# end if 2 if ( not $glob_reached_optimal_h) then # if number 2 $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 2 hnew = sz2 #END block return(hnew) #BOTTOM ADJUST FOR POLE end # End Function number 9 # Begin Function number 10 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 2 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 2 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 10 # Begin Function number 11 def check_for_pole() #TOP CHECK FOR POLE $array_pole[1] = $glob_large_float $array_pole[2] = $glob_large_float tmp_rad = $glob_large_float prev_tmp_rad = $glob_large_float tmp_ratio = $glob_large_float rad_c = $glob_large_float $array_poles[1][1] = $glob_large_float $array_poles[1][2] = $glob_large_float #TOP radius ratio test in Henrici1 found_sing = 1 n = $glob_max_terms - 1 - 10 cnt = 0 while ((cnt < 5) and (found_sing == 1)) do # do number 1 if ((omniabs($array_y_higher[1][n]) == 0.0) or (omniabs($array_y_higher[1][n+1]) == 0.0)) then # if number 2 found_sing = 0 else tmp_rad = omniabs($array_y_higher[1][n] * $glob_h / $array_y_higher[1][n + 1]) tmp_ratio = tmp_rad / prev_tmp_rad if ((cnt > 0 ) and (tmp_ratio < 2.0) and (tmp_ratio > 0.5)) then # if number 3 if (tmp_rad < rad_c) then # if number 4 rad_c = tmp_rad end# end if 4 elsif (cnt == 0) then # if number 4 if (tmp_rad < rad_c) then # if number 5 rad_c = tmp_rad end# end if 5 elsif (cnt > 0) then # if number 5 found_sing = 0 end# end if 5 end# end if 4 prev_tmp_rad = tmp_rad cnt = cnt + 1 n = n + 1 end# end do number 1 if (found_sing == 1) then # if number 4 if (rad_c < $array_pole[1]) then # if number 5 $array_pole[1] = rad_c $array_poles[1][1] = rad_c end# end if 5 end# end if 4 #BOTTOM radius ratio test in Henrici1 #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]) == 0.0) or (omniabs($array_y_higher[1][m-1]) == 0.0) or (omniabs($array_y_higher[1][m-2]) == 0.0))) do # do number 1 m = m - 1 end# end do number 1 if (m > 10) then # if number 4 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) > 0.0) then # if number 5 rcs = $glob_h/hdrc ord_no = (rm1*convfloat((m-2)*(m-2))-rm0*convfloat(m-3))/hdrc $array_real_poles[1][1] = rcs $array_real_poles[1][2] = ord_no else $array_real_poles[1][1] = $glob_large_float $array_real_poles[1][2] = $glob_large_float end# end if 5 else $array_real_poles[1][1] = $glob_large_float $array_real_poles[1][2] = $glob_large_float end# end if 4 #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 1 if (omniabs($array_y_higher[1][n]) != 0.0) then # if number 4 cnt = cnt + 1 else cnt = 0 end# end if 4 n = n - 1 end# end do number 1 m = n + cnt if (m <= 10) then # if number 4 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) == 0.0) or (omniabs(dr1) == 0.0)) then # if number 5 rad_c = $glob_large_float ord_no = $glob_large_float else if (omniabs(nr1*dr2 - nr2 * dr1) != 0.0) then # if number 6 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) != 0.0) then # if number 7 if (rcs > 0.0) then # if number 8 rad_c = sqrt(rcs) * omniabs($glob_h) else rad_c = $glob_large_float end# end if 8 else rad_c = $glob_large_float ord_no = $glob_large_float end# end if 7 else rad_c = $glob_large_float ord_no = $glob_large_float end# end if 6 end# end if 5 $array_complex_poles[1][1] = rad_c $array_complex_poles[1][2] = ord_no end# end if 4 #BOTTOM RADII COMPLEX EQ = 1 #START ADJUST ALL SERIES if ($array_pole[1] * $glob_ratio_of_radius < omniabs($glob_h)) then # if number 4 h_new = $array_pole[1] * $glob_ratio_of_radius term = 1 ratio = 1.0 while (term <= $glob_max_terms) do # do number 1 $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 1 $glob_h = h_new end# end if 4 #BOTTOM ADJUST ALL SERIES if (reached_interval()) then # if number 4 display_poles() end# end if 4 end # End Function number 11 # Begin Function number 12 def get_norms() if ( not $glob_initial_pass) then # if number 4 iii = 1 while (iii <= $glob_max_terms) do # do number 1 $array_norms[iii] = 0.0 iii = iii + 1 end# end do number 1 #TOP GET NORMS iii = 1 while (iii <= $glob_max_terms) do # do number 1 if (omniabs($array_y[iii]) > $array_norms[iii]) then # if number 5 $array_norms[iii] = omniabs($array_y[iii]) end# end if 5 iii = iii + 1 end# end do number 1 #BOTTOM GET NORMS end# end if 4 end # End Function number 12 # Begin Function number 13 def atomall() #TOP ATOMALL #END OUTFILE1 #BEGIN ATOMHDR1 #emit pre mult CONST - LINEAR $eq_no = 1 i = 1 $array_tmp1[1] = $array_const_0D3[1] * $array_x[1] #emit pre add LINEAR - CONST $eq_no = 1 i = 1 $array_tmp2[1] = $array_tmp1[1] + $array_const_0D1[1] #emit pre sin 1 $eq_no = 1 $array_tmp3[1] = sin($array_tmp2[1]) $array_tmp3_g[1] = cos($array_tmp2[1]) #emit pre add CONST FULL $eq_no = 1 i = 1 $array_tmp4[1] = $array_const_0D0[1] + $array_tmp3[1] #emit pre mult CONST - LINEAR $eq_no = 1 i = 1 $array_tmp5[1] = $array_const_0D1[1] * $array_x[1] #emit pre add LINEAR - CONST $eq_no = 1 i = 1 $array_tmp6[1] = $array_tmp5[1] + $array_const_0D2[1] #emit pre add FULL - LINEAR $eq_no = 1 i = 1 $array_tmp7[1] = $array_tmp4[1] + $array_tmp6[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_tmp7[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_0D3[1] * $array_x[2] #emit pre add LINEAR - CONST $eq_no = 1 i = 2 $array_tmp2[2] = $array_tmp1[2] #emit pre sin ID_LINEAR iii = 2 $eq_no = 1 $array_tmp3[2] = $array_tmp3_g[1] * $array_tmp2[2] / 1 $array_tmp3_g[2] = -$array_tmp3[1] * $array_tmp2[2] / 1 #emit pre add CONST FULL $eq_no = 1 i = 2 $array_tmp4[2] = $array_tmp3[2] #emit pre mult CONST - LINEAR $eq_no = 1 i = 2 $array_tmp5[2] = $array_const_0D1[1] * $array_x[2] #emit pre add LINEAR - CONST $eq_no = 1 i = 2 $array_tmp6[2] = $array_tmp5[2] #emit pre add FULL - LINEAR $eq_no = 1 i = 2 $array_tmp7[2] = $array_tmp4[2] + $array_tmp6[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_tmp7[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 sin ID_LINEAR iii = 3 $eq_no = 1 $array_tmp3[3] = $array_tmp3_g[2] * $array_tmp2[2] / 2 $array_tmp3_g[3] = -$array_tmp3[2] * $array_tmp2[2] / 2 #emit pre add CONST FULL $eq_no = 1 i = 3 $array_tmp4[3] = $array_tmp3[3] #emit pre add FULL - LINEAR $eq_no = 1 i = 3 $array_tmp7[3] = $array_tmp4[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_tmp7[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 sin ID_LINEAR iii = 4 $eq_no = 1 $array_tmp3[4] = $array_tmp3_g[3] * $array_tmp2[2] / 3 $array_tmp3_g[4] = -$array_tmp3[3] * $array_tmp2[2] / 3 #emit pre add CONST FULL $eq_no = 1 i = 4 $array_tmp4[4] = $array_tmp3[4] #emit pre add FULL - LINEAR $eq_no = 1 i = 4 $array_tmp7[4] = $array_tmp4[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_tmp7[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 sin ID_LINEAR iii = 5 $eq_no = 1 $array_tmp3[5] = $array_tmp3_g[4] * $array_tmp2[2] / 4 $array_tmp3_g[5] = -$array_tmp3[4] * $array_tmp2[2] / 4 #emit pre add CONST FULL $eq_no = 1 i = 5 $array_tmp4[5] = $array_tmp3[5] #emit pre add FULL - LINEAR $eq_no = 1 i = 5 $array_tmp7[5] = $array_tmp4[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_tmp7[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 sin LINEAR $eq_no = 1 $array_tmp3[kkk] = $array_tmp3_g[kkk - 1] * $array_tmp2[2] / (kkk - 1) $array_tmp3_g[kkk] = -$array_tmp3[kkk - 1] * $array_tmp2[2] / (kkk - 1) #emit NOT FULL - FULL add $eq_no = 1 $array_tmp4[kkk] = $array_tmp3[kkk] #emit FULL - NOT FULL add $eq_no = 1 $array_tmp7[kkk] = $array_tmp4[kkk] #emit assign $eq_no = 1 order_d = 1 if (kkk + order_d < $glob_max_terms) then # if number 1 if ( not $array_y_set_initial[1][kkk + order_d]) then # if number 2 temporary = $array_tmp7[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 1 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 1 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 13 #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,m,radius,order2) if (typ == 1) then # if number 8 omniout_str(ALWAYS,"Real") else omniout_str(ALWAYS,"Complex") end# end if 8 omniout_int(ALWAYS,"m",4, m ,4," ") 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 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") elsif pole == 4 then # if number 11 file.printf("Yes") else file.printf("No") end# end if 11 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 11 omniout_str(ALWAYS,"Illegal max_terms = -- Using 30") $glob_max_terms = 30 end# end if 11 if ($glob_max_iter < 2) then # if number 11 omniout_str(ALWAYS,"Illegal max_iter") errflag = true end# end if 11 if (errflag) then # if number 11 end# end if 11 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 11 sec_left = 0.0 else if (sub2 > 0.0) then # if number 12 rrr = (sub1/sub2) sec_left = rrr * ms2 - ms2 else sec_left = 0.0 end# end if 12 end# end if 11 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 11 rrr = (100.0*sub2)/sub1 else rrr = 0.0 end# end if 11 return ( rrr) end # End Function number 27 # Begin Function number 28 def factorial_2(nnn) if (nnn > 0) then # if number 11 return ( nnn * factorial_2(nnn - 1)) else return ( 1.0) end# end if 11 end # End Function number 28 # Begin Function number 29 def factorial_1(nnn) if (nnn <= $glob_max_terms) then # if number 11 if ($array_fact_1[nnn] == 0) then # if number 12 ret = factorial_2(nnn) $array_fact_1[nnn] = ret else ret = $array_fact_1[nnn] end# end if 12 else ret = factorial_2(nnn) end# end if 11 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 11 if ($array_fact_2[mmm][nnn] == 0) then # if number 12 ret = factorial_1(mmm)/factorial_1(nnn) $array_fact_2[mmm][nnn] = ret else ret = $array_fact_2[mmm][nnn] end# end if 12 else ret = factorial_2(mmm)/factorial_2(nnn) end# end if 11 return ( ret) end # End Function number 30 # Begin Function number 31 def convfp(mmm) return(mmm.to_f) end # End Function number 31 # Begin Function number 32 def convfloat(mmm) return(mmm.to_f) 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) if ((x <= 0.0) and (y < 0.0)) then # if number 11 puts "expt error x = " + x.to_s + "y = " + y.to_s end# end if 11 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(0.05 * x * x + 0.2 * x - cos(0.3 * x + 0.1) / 0.3) 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_estimated_step_error = 0.0 $glob_ratio_of_radius = 0.1 $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_min_h = 0.000001 $glob_type_given_pole = 0 $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.0 $glob_smallish_float = 0.0 $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/add_full_linpostode.ode#################") omniout_str(ALWAYS,"diff ( y , x , 1 ) = sin(0.3 * x + 0.1) + (0.1 * x + 0.2) ;") 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_display_interval=0.1 ") omniout_str(ALWAYS,"$glob_look_poles=true ") omniout_str(ALWAYS,"$glob_max_iter=1000000 ") 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.01 ") 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(0.05 * x * x + 0.2 * x - cos(0.3 * x + 0.1) / 0.3) ") 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 = 0.0 $glob_smallish_float = 0.0 $glob_large_float = 1.0e100 $glob_almost_1 = 0.99 #BEGIN FIRST INPUT BLOCK #BEGIN FIRST INPUT BLOCK # Digits:=32; ELIMINATED in preodein.rb max_terms=30 #END FIRST INPUT BLOCK #START OF INITS AFTER INPUT BLOCK $glob_max_terms = max_terms $glob_html_log = true #END OF INITS AFTER INPUT BLOCK $array_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(4 + 1) $array_real_pole= Array.new(4 + 1) $array_complex_pole= Array.new(4 + 1) $array_1st_rel_error= Array.new(2 + 1) $array_last_rel_error= Array.new(2 + 1) $array_type_pole= Array.new(2 + 1) $array_type_real_pole= Array.new(2 + 1) $array_type_complex_pole= Array.new(2 + 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_g= 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_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(2 + 1 ,3 + 1 ) $array_given_rad_poles = array2d(2 + 1 ,3 + 1 ) $array_given_ord_poles = array2d(2 + 1 ,3 + 1 ) $array_real_poles = array2d(2 + 1 ,3 + 1 ) $array_complex_poles = array2d(2 + 1 ,3 + 1 ) $array_fact_2 = array2d(max_terms + 1 ,max_terms + 1 ) term = 1 while (term <= max_terms) do # do number 1 $array_y_init[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_norms[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_fact_1[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 4) do # do number 1 $array_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 4) do # do number 1 $array_real_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 4) do # do number 1 $array_complex_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 2) do # do number 1 $array_1st_rel_error[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 2) do # do number 1 $array_last_rel_error[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 2) do # do number 1 $array_type_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 2) do # do number 1 $array_type_real_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= 2) do # do number 1 $array_type_complex_pole[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_y[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_x[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp0[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp1[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp2[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp3_g[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp3[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp4[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp5[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp6[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_tmp7[term] = 0.0 term = term + 1 end# end do number 1 term = 1 while (term <= max_terms) do # do number 1 $array_m1[term] = 0.0 term = term + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= max_terms) do # do number 2 $array_y_higher[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= max_terms) do # do number 2 $array_y_higher_work[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= max_terms) do # do number 2 $array_y_higher_work2[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= max_terms) do # do number 2 $array_y_set_initial[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= 3) do # do number 2 $array_poles[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= 3) do # do number 2 $array_given_rad_poles[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= 3) do # do number 2 $array_given_ord_poles[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= 3) do # do number 2 $array_real_poles[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=2) do # do number 1 term = 1 while (term <= 3) do # do number 2 $array_complex_poles[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 ord = 1 while (ord <=max_terms) do # do number 1 term = 1 while (term <= max_terms) do # do number 2 $array_fact_2[ord][term] = 0.0 term = term + 1 end# end do number 2 ord = ord + 1 end# end do number 1 #BEGIN ARRAYS DEFINED AND INITIALIZATED $array_y = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_y[term] = 0.0 term = term + 1 end# end do number 1 $array_x = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_x[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp0 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp0[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp1[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp2 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp2[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp3_g = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp3_g[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp3 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp3[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp4 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp4[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp5 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp5[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp6 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp6[term] = 0.0 term = term + 1 end# end do number 1 $array_tmp7 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_tmp7[term] = 0.0 term = term + 1 end# end do number 1 $array_m1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_m1[term] = 0.0 term = term + 1 end# end do number 1 $array_const_1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_const_1[term] = 0.0 term = term + 1 end# end do number 1 $array_const_1[1] = 1 $array_const_0D0 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_const_0D0[term] = 0.0 term = term + 1 end# end do number 1 $array_const_0D0[1] = 0.0 $array_const_0D3 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_const_0D3[term] = 0.0 term = term + 1 end# end do number 1 $array_const_0D3[1] = 0.3 $array_const_0D1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms + 1) do # do number 1 $array_const_0D1[term] = 0.0 term = term + 1 end# end do number 1 $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 1 $array_const_0D2[term] = 0.0 term = term + 1 end# end do number 1 $array_const_0D2[1] = 0.2 $array_m1 = Array.new(max_terms+1 + 1) term = 1 while (term <= max_terms) do # do number 1 $array_m1[term] = 0.0 term = term + 1 end# end do number 1 $array_m1[1] = -1.0 #END ARRAYS DEFINED AND INITIALIZATED #Initing Factorial Tables iiif = 0 while (iiif <= $glob_max_terms) do # do number 1 jjjf = 0 while (jjjf <= $glob_max_terms) do # do number 2 $array_fact_1[iiif] = 0 $array_fact_2[iiif][jjjf] = 0 jjjf = jjjf + 1 end# end do number 2 iiif = iiif + 1 end# end do number 1 #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_display_interval=0.1 $glob_look_poles=true $glob_max_iter=1000000 $glob_max_minutes=10 #END SECOND INPUT BLOCK #BEGIN OVERRIDE BLOCK $glob_desired_digits_correct=10 $glob_display_interval=0.01 $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 #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) found_h = false $glob_h = $glob_min_h if ($glob_max_h < $glob_h) then # if number 4 $glob_h = $glob_max_h end# end if 4 if ($glob_display_interval < $glob_h) then # if number 4 $glob_h = $glob_display_interval end# end if 4 best_h = $glob_h min_value = $glob_large_float est_answer = est_size_answer() opt_iter = 1 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,"") estimated_step_error = 0.0 while ((opt_iter <= 100) and ( not found_h)) do # do number 1 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 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 atomall() estimated_step_error = test_suggested_h() omniout_float(ALWAYS,"estimated_step_error",32,estimated_step_error,32,"") if (((estimated_step_error > est_needed_step_err) and (opt_iter == 1)) or ($glob_h >= $glob_max_h )) then # if number 4 found_h = true $glob_h = $glob_max_h best_h = $glob_h elsif ((estimated_step_error > est_needed_step_err) and ( not found_h)) then # if number 5 $glob_h = $glob_h/2.0 best_h = $glob_h found_h = true else $glob_h = $glob_h*2.0 best_h = $glob_h end# end if 5 omniout_float(ALWAYS,"best_h",32,best_h,32,"") opt_iter = opt_iter + 1 end# end do number 1 if (( not found_h) and (opt_iter == 1)) then # if number 5 omniout_str(ALWAYS,"Beginning $glob_h too large.") found_h = false end# end if 5 if (opt_iter > 100) then # if number 5 $glob_h = $glob_max_h found_h = false end# end if 5 if ($glob_display_interval < $glob_h) then # if number 5 $glob_h = $glob_display_interval end# end if 5 #END OPTIMIZE CODE if ($glob_html_log) then # if number 5 html_log_file = File.new("entry.html","w") end# end if 5 #BEGIN SOLUTION CODE if (found_h) then # if number 5 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 1 $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 1 rows = order_diff r_order = 1 while (r_order <= rows) do # do number 1 term_no = 1 while (term_no <= (rows - r_order + 1)) do # do number 2 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 2 r_order = r_order + 1 end# end do number 1 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 1 #left paren 0001C if (reached_interval()) then # if number 6 omniout_str(INFO," ") omniout_str(INFO,"TOP MAIN SOLVE Loop") end# end if 6 $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 6 #left paren 0004C check_for_pole() end# end if 6#was right paren 0004C if (reached_interval()) then # if number 6 $glob_next_display = $glob_next_display + $glob_display_interval end# end if 6 $array_x[1] = $array_x[1] + $glob_h $array_x[2] = $glob_h #Jump Series $array_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 2 $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 2 #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 2 temp_sum = temp_sum + $array_y_higher_work[ord][iii] iii = iii - 1 end# end do number 2 $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 2 $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 2 #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 2 temp_sum = temp_sum + $array_y_higher_work[ord][iii] iii = iii - 1 end# end do number 2 $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 2 $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 2 #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 2 temp_sum = temp_sum + $array_y_higher_work[ord][iii] iii = iii - 1 end# end do number 2 $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 2 $array_y[term_no] = $array_y_higher_work2[1][term_no] ord = 1 while (ord <= order_diff) do # do number 3 $array_y_higher[ord][term_no] = $array_y_higher_work2[ord][term_no] ord = ord + 1 end# end do number 3 term_no = term_no - 1 end# end do number 2 #END PART 2 HEVE MOVED TERMS to REGULAR Array end# end do number 1#right paren 0001C omniout_str(ALWAYS,"Finished!") if ($glob_iter >= $glob_max_iter) then # if number 6 omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!") end# end if 6 if (elapsed_time_seconds() - convfloat($glob_orig_start_sec) >= convfloat($glob_max_sec )) then # if number 6 omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!") end# end if 6 $glob_clock_sec = elapsed_time_seconds() omniout_str(INFO,"diff ( y , x , 1 ) = sin(0.3 * x + 0.1) + (0.1 * x + 0.2) ;") omniout_int(INFO,"Iterations ",32,$glob_iter,4," ") prog_report(x_start,x_end) if ($glob_html_log) then # if number 6 logstart(html_log_file) logitem_str(html_log_file,"2013-05-25T23:41:45-05:00") logitem_str(html_log_file,"Ruby") logitem_str(html_log_file,"add_full_lin") logitem_str(html_log_file,"diff ( y , x , 1 ) = sin(0.3 * x + 0.1) + (0.1 * x + 0.2) ;") 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_time(html_log_file,convfloat($glob_clock_sec)) if ($glob_percent_done < 100.0) then # if number 7 logitem_time(html_log_file,convfloat($glob_total_exp_sec)) 0 else logitem_str(html_log_file,"Done") 0 end# end if 7 log_revs(html_log_file," 189 ") logitem_str(html_log_file,"add_full_lin diffeq.rb") logitem_str(html_log_file,"add_full_lin Ruby results") logitem_str(html_log_file,"All Tests - All Languages") logend(html_log_file) end# end if 6 if ($glob_html_log) then # if number 6 html_log_file.close; end# end if 6 end# end if 5 #END OUTFILEMAIN end # End Function number 13 main()