|\^/| Maple 12 (IBM INTEL LINUX) ._|\| |/|_. Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2008 \ MAPLE / All rights reserved. Maple is a trademark of <____ ____> Waterloo Maple Inc. | Type ? for help. > #BEGIN OUTFILE1 > > # Begin Function number 3 > reached_interval := proc() > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > > > local ret; > > if ((((array_x[1] >= glob_next_display) and not glob_neg_h) or ((array_x[1] <= glob_next_display) and glob_neg_h)) or (glob_next_display = 0.0)) then # if number 1 > ret := true; > else > ret := false; > fi;# end if 1 > ; > return(ret); > > # End Function number 3 > end; reached_interval := proc() local ret; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; if glob_next_display <= array_x[1] and not glob_neg_h or array_x[1] <= glob_next_display and glob_neg_h or glob_next_display = 0. then ret := true else ret := false end if; return ret end proc > # Begin Function number 4 > display_alot := proc(iter) > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local abserr, analytic_val_y, ind_var, numeric_val, relerr, term_no; > > > > > > #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.0) then # if number 4 > glob_good_digits := -trunc(log10(relerr/100.0)); > else > glob_good_digits := Digits; > fi;# end if 4 > ; > else > relerr := -1.0 ; > glob_good_digits := -1; > fi;# end if 3 > ; > if (glob_iter = 1) then # if number 3 > array_1st_rel_error[1] := relerr; > else > array_last_rel_error[1] := relerr; > fi;# 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," "); > fi;# end if 2 > ; > #BOTTOM DISPLAY ALOT > fi;# end if 1 > ; > > # End Function number 4 > end; display_alot := proc(iter) local abserr, analytic_val_y, ind_var, numeric_val, relerr, term_no; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; if reached_interval() then if 0 <= iter then 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. then relerr := abserr*100.0/omniabs(analytic_val_y); if relerr <> 0. then glob_good_digits := -trunc(log10(relerr/100.0)) else glob_good_digits := Digits end if else relerr := -1.0; glob_good_digits := -1 end if; if glob_iter = 1 then array_1st_rel_error[1] := relerr else array_last_rel_error[1] := relerr end if; 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 if end if end proc > # Begin Function number 5 > adjust_for_pole := proc(h_param) > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local hnew, sz2, tmp; > > > > #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; > fi;# end if 2 > fi;# 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); > fi;# end if 2 > fi;# 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]; > fi;# end if 1 > ; > hnew := sz2; > ;#END block > return(hnew); > #BOTTOM ADJUST FOR POLE > > # End Function number 5 > end; adjust_for_pole := proc(h_param) local hnew, sz2, tmp; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; hnew := h_param; glob_normmax := glob_small_float; if glob_small_float < omniabs(array_y_higher[1, 1]) then tmp := omniabs(array_y_higher[1, 1]); if tmp < glob_normmax then glob_normmax := tmp end if end if; if glob_look_poles and glob_small_float < omniabs(array_pole[1]) and array_pole[1] <> glob_large_float then sz2 := array_pole[1]/10.0; if sz2 < hnew then omniout_float(INFO, "glob_h adjusted to ", 20, h_param, 12, "due to singularity."); omniout_str(INFO, "Reached Optimal"); return hnew end if end if; if not glob_reached_optimal_h then 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 if; hnew := sz2; return hnew end proc > # Begin Function number 6 > prog_report := proc(x_start,x_end) > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local clock_sec, opt_clock_sec, clock_sec1, expect_sec, left_sec, percent_done, total_clock_sec; > > > > > > #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)); > 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)); > fi;# 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 Function number 6 > end; prog_report := proc(x_start, x_end) local clock_sec, opt_clock_sec, clock_sec1, expect_sec, left_sec, percent_done, total_clock_sec; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; 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)); 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 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)) end if; omniout_str_noeol(INFO, "Time to Timeout "); omniout_timestr(convfloat(left_sec)); omniout_float(INFO, "Percent Done ", 33, percent_done, 4, "%") end proc > # Begin Function number 7 > check_for_pole := proc() > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local cnt, dr1, dr2, ds1, ds2, hdrc, m, n, nr1, nr2, ord_no, rad_c, rcs, rm0, rm1, rm2, rm3, rm4, found; > > > > > > #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; > od;# 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; > fi;# end if 2 > else > array_real_pole[1,1] := glob_large_float; > array_real_pole[1,2] := glob_large_float; > fi;# 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; > fi;# end if 1 > ; > n := n - 1; > od;# end do number 2 > ; > m := n + cnt; > if (m <= 10) then # if number 1 > array_complex_pole[1,1] := glob_large_float; > array_complex_pole[1,2] := glob_large_float; > elif ((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))) then # if number 2 > array_complex_pole[1,1] := glob_large_float; > array_complex_pole[1,2] := 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 > array_complex_pole[1,1] := glob_large_float; > array_complex_pole[1,2] := 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; > fi;# end if 6 > else > rad_c := glob_large_float; > ord_no := glob_large_float; > fi;# end if 5 > else > rad_c := glob_large_float; > ord_no := glob_large_float; > fi;# end if 4 > fi;# end if 3 > ; > array_complex_pole[1,1] := rad_c; > array_complex_pole[1,2] := ord_no; > fi;# 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 > omniout_str(ALWAYS,"Complex estimate of poles used"); > fi;# end if 3 > ; > fi;# 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 > omniout_str(ALWAYS,"Real estimate of pole used"); > fi;# end if 3 > ; > fi;# 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"); > fi;# end if 3 > ; > fi;# 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 > omniout_str(ALWAYS,"Real estimate of pole used"); > fi;# end if 3 > ; > fi;# 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 > omniout_str(ALWAYS,"Complex estimate of poles used"); > fi;# end if 3 > ; > fi;# 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"); > fi;# end if 3 > ; > fi;# 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]; > fi;# end if 2 > ; > #BOTTOM WHICH RADIUS EQ = 1 > #BOTTOM CHECK FOR POLE > if (reached_interval()) then # if number 2 > display_pole(); > fi;# end if 2 > > # End Function number 7 > end; check_for_pole := proc() local cnt, dr1, dr2, ds1, ds2, hdrc, m, n, nr1, nr2, ord_no, rad_c, rcs, rm0, rm1, rm2, rm3, rm4, found; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; n := glob_max_terms; m := n - 2; while 10 <= m 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 m := m - 1 end do; if 10 < m then 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 glob_small_float < omniabs(hdrc) then 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 if else array_real_pole[1, 1] := glob_large_float; array_real_pole[1, 2] := glob_large_float end if; n := glob_max_terms - 2; cnt := 0; while cnt < 5 and 10 <= n do if glob_small_float < omniabs(array_y_higher[1, n]) then cnt := cnt + 1 else cnt := 0 end if; n := n - 1 end do; m := n + cnt; if m <= 10 then array_complex_pole[1, 1] := glob_large_float; array_complex_pole[1, 2] := glob_large_float elif glob_large_float <= omniabs(array_y_higher[1, m]) or glob_large_float <= omniabs(array_y_higher[1, m - 1]) or glob_large_float <= omniabs(array_y_higher[1, m - 2]) or glob_large_float <= omniabs(array_y_higher[1, m - 3]) or glob_large_float <= omniabs(array_y_higher[1, m - 4]) or glob_large_float <= omniabs(array_y_higher[1, m - 5]) then array_complex_pole[1, 1] := glob_large_float; array_complex_pole[1, 2] := 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)*(1.0)/rm1 + 2.0/rm2 - 1.0/rm3; dr2 := (-1)*(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 array_complex_pole[1, 1] := glob_large_float; array_complex_pole[1, 2] := glob_large_float else if glob_small_float < omniabs(nr1*dr2 - nr2*dr1) then rcs := (ds1*dr2 - ds2*dr1 + dr1*dr2)/(nr1*dr2 - nr2*dr1); ord_no := (rcs*nr1 - ds1)/(2.0*dr1) - convfloat(m)/2.0; if glob_small_float < omniabs(rcs) then if 0. < rcs then rad_c := sqrt(rcs)*omniabs(glob_h) else rad_c := glob_large_float end if else rad_c := glob_large_float; ord_no := glob_large_float end if else rad_c := glob_large_float; ord_no := glob_large_float end if end if; array_complex_pole[1, 1] := rad_c; array_complex_pole[1, 2] := ord_no end if; found := false; 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 0. < array_complex_pole[1, 1] and 0. < array_complex_pole[1, 2] then 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 omniout_str(ALWAYS, "Complex estimate of poles used") end if end if; if not found and array_real_pole[1, 1] <> glob_large_float and array_real_pole[1, 2] <> glob_large_float and 0. < array_real_pole[1, 1] and 0. < array_real_pole[1, 2] 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. or array_complex_pole[1, 2] <= 0.) then 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 omniout_str(ALWAYS, "Real estimate of pole used") end if end if; 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 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 omniout_str(ALWAYS, "NO POLE") end if end if; if not found and array_real_pole[1, 1] < array_complex_pole[1, 1] and 0. < array_real_pole[1, 1] and 0. < array_real_pole[1, 2] then 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 omniout_str(ALWAYS, "Real estimate of pole used") end if end if; if not found and array_complex_pole[1, 1] <> glob_large_float and array_complex_pole[1, 2] <> glob_large_float and 0. < array_complex_pole[1, 1] and 0. < array_complex_pole[1, 2] then 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 omniout_str(ALWAYS, "Complex estimate of poles used") end if end if; if not found then array_poles[1, 1] := glob_large_float; array_poles[1, 2] := glob_large_float; array_type_pole[1] := 3; if reached_interval() then omniout_str(ALWAYS, "NO POLE") end if end if; array_pole[1] := glob_large_float; array_pole[2] := glob_large_float; if array_poles[1, 1] < array_pole[1] then array_pole[1] := array_poles[1, 1]; array_pole[2] := array_poles[1, 2] end if; if reached_interval() then display_pole() end if end proc > # Begin Function number 8 > get_norms := proc() > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local iii; > > > > 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; > od;# 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]); > fi;# end if 3 > ; > iii := iii + 1; > od;# end do number 2 > #BOTTOM GET NORMS > ; > fi;# end if 2 > ; > > # End Function number 8 > end; get_norms := proc() local iii; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; if not glob_initial_pass then iii := 1; while iii <= glob_max_terms do array_norms[iii] := 0.; iii := iii + 1 end do; iii := 1; while iii <= glob_max_terms do if array_norms[iii] < omniabs(array_y[iii]) then array_norms[iii] := omniabs(array_y[iii]) end if; iii := iii + 1 end do end if end proc > # Begin Function number 9 > atomall := proc() > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > > local kkk, order_d, adj2, temporary, term; > > > > > > #TOP ATOMALL > #END OUTFILE1 > #BEGIN ATOMHDR1 > #emit pre mult FULL CONST $eq_no = 1 i = 1 > array_tmp1[1] := array_m1[1] * array_const_2D0[1]; > #emit pre mult FULL LINEAR $eq_no = 1 i = 1 > #emit pre mult LINEAR - FULL $eq_no = 1 i = 1 > array_tmp2[1] := array_tmp1[1] * array_x[1]; > #emit pre mult LINEAR - LINEAR $eq_no = 1 i = 1 > array_tmp3[1] := array_x[1] * array_x[1]; > #emit pre add FULL - CONST $eq_no = 1 i = 1 > array_tmp4[1] := array_tmp3[1] + array_const_1D0[1]; > #emit pre div FULL - FULL $eq_no = 1 i = 1 > array_tmp5[1] := (array_tmp2[1] / (array_tmp4[1])); > #emit pre mult LINEAR - LINEAR $eq_no = 1 i = 1 > array_tmp6[1] := array_x[1] * array_x[1]; > #emit pre add FULL - CONST $eq_no = 1 i = 1 > array_tmp7[1] := array_tmp6[1] + array_const_1D0[1]; > #emit pre div FULL - FULL $eq_no = 1 i = 1 > array_tmp8[1] := (array_tmp5[1] / (array_tmp7[1])); > #emit pre add CONST FULL $eq_no = 1 i = 1 > array_tmp9[1] := array_const_0D0[1] + array_tmp8[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_tmp9[1] * expt(glob_h , (1)) * factorial_3(0,1); > array_y[2] := temporary; > array_y_higher[1,2] := temporary; > temporary := temporary / glob_h * (2.0); > array_y_higher[2,1] := temporary > ; > fi;# end if 2 > ; > fi;# end if 1 > ; > kkk := 2; > #END ATOMHDR1 > #BEGIN ATOMHDR2 > #emit pre mult FULL CONST $eq_no = 1 i = 2 > array_tmp1[2] := array_m1[2] * array_const_2D0[1]; > #emit pre mult LINEAR FULL $eq_no = 1 i = 2 > array_tmp2[2] := array_tmp1[2] * array_x[kkk - 1] + array_tmp1[1] * array_x[kkk]; > #emit pre mult LINEAR - LINEAR $eq_no = 1 i = 2 > array_tmp3[2] := array_x[1] * array_x[2] + array_x[2] * array_x[1]; > #emit pre add FULL CONST $eq_no = 1 i = 2 > array_tmp4[2] := array_tmp3[2]; > #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 mult LINEAR - LINEAR $eq_no = 1 i = 2 > array_tmp6[2] := array_x[1] * array_x[2] + array_x[2] * array_x[1]; > #emit pre add FULL CONST $eq_no = 1 i = 2 > array_tmp7[2] := array_tmp6[2]; > #emit pre div FULL - FULL $eq_no = 1 i = 2 > array_tmp8[2] := ((array_tmp5[2] - ats(2,array_tmp7,array_tmp8,2))/array_tmp7[1]); > #emit pre add CONST FULL $eq_no = 1 i = 2 > array_tmp9[2] := array_tmp8[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_tmp9[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 > ; > fi;# end if 2 > ; > fi;# end if 1 > ; > kkk := 3; > #END ATOMHDR2 > #BEGIN ATOMHDR3 > #emit pre mult FULL CONST $eq_no = 1 i = 3 > array_tmp1[3] := array_m1[3] * array_const_2D0[1]; > #emit pre mult LINEAR FULL $eq_no = 1 i = 3 > array_tmp2[3] := array_tmp1[2] * array_x[kkk - 1] + array_tmp1[1] * array_x[kkk]; > #emit pre mult LINEAR - LINEAR $eq_no = 1 i = 3 > array_tmp3[3] := array_x[2] * array_x[2]; > #emit pre add FULL CONST $eq_no = 1 i = 3 > array_tmp4[3] := array_tmp3[3]; > #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 mult LINEAR - LINEAR $eq_no = 1 i = 3 > array_tmp6[3] := array_x[2] * array_x[2]; > #emit pre add FULL CONST $eq_no = 1 i = 3 > array_tmp7[3] := array_tmp6[3]; > #emit pre div FULL - FULL $eq_no = 1 i = 3 > array_tmp8[3] := ((array_tmp5[3] - ats(3,array_tmp7,array_tmp8,2))/array_tmp7[1]); > #emit pre add CONST FULL $eq_no = 1 i = 3 > array_tmp9[3] := array_tmp8[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_tmp9[3] * expt(glob_h , (1)) * factorial_3(2,3); > array_y[4] := temporary; > array_y_higher[1,4] := temporary; > temporary := temporary / glob_h * (2.0); > array_y_higher[2,3] := temporary > ; > fi;# end if 2 > ; > fi;# end if 1 > ; > kkk := 4; > #END ATOMHDR3 > #BEGIN ATOMHDR4 > #emit pre mult FULL CONST $eq_no = 1 i = 4 > array_tmp1[4] := array_m1[4] * array_const_2D0[1]; > #emit pre mult LINEAR FULL $eq_no = 1 i = 4 > array_tmp2[4] := array_tmp1[2] * array_x[kkk - 1] + array_tmp1[1] * array_x[kkk]; > #emit pre add FULL CONST $eq_no = 1 i = 4 > array_tmp4[4] := array_tmp3[4]; > #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 FULL CONST $eq_no = 1 i = 4 > array_tmp7[4] := array_tmp6[4]; > #emit pre div FULL - FULL $eq_no = 1 i = 4 > array_tmp8[4] := ((array_tmp5[4] - ats(4,array_tmp7,array_tmp8,2))/array_tmp7[1]); > #emit pre add CONST FULL $eq_no = 1 i = 4 > array_tmp9[4] := array_tmp8[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_tmp9[4] * expt(glob_h , (1)) * factorial_3(3,4); > array_y[5] := temporary; > array_y_higher[1,5] := temporary; > temporary := temporary / glob_h * (2.0); > array_y_higher[2,4] := temporary > ; > fi;# end if 2 > ; > fi;# end if 1 > ; > kkk := 5; > #END ATOMHDR4 > #BEGIN ATOMHDR5 > #emit pre mult FULL CONST $eq_no = 1 i = 5 > array_tmp1[5] := array_m1[5] * array_const_2D0[1]; > #emit pre mult LINEAR FULL $eq_no = 1 i = 5 > array_tmp2[5] := array_tmp1[2] * array_x[kkk - 1] + array_tmp1[1] * array_x[kkk]; > #emit pre add FULL CONST $eq_no = 1 i = 5 > array_tmp4[5] := array_tmp3[5]; > #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 FULL CONST $eq_no = 1 i = 5 > array_tmp7[5] := array_tmp6[5]; > #emit pre div FULL - FULL $eq_no = 1 i = 5 > array_tmp8[5] := ((array_tmp5[5] - ats(5,array_tmp7,array_tmp8,2))/array_tmp7[1]); > #emit pre add CONST FULL $eq_no = 1 i = 5 > array_tmp9[5] := array_tmp8[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_tmp9[5] * expt(glob_h , (1)) * factorial_3(4,5); > array_y[6] := temporary; > array_y_higher[1,6] := temporary; > temporary := temporary / glob_h * (2.0); > array_y_higher[2,5] := temporary > ; > fi;# end if 2 > ; > fi;# end if 1 > ; > kkk := 6; > #END ATOMHDR5 > #BEGIN OUTFILE3 > #Top Atomall While Loop-- outfile3 > while (kkk <= glob_max_terms) do # do number 1 > #END OUTFILE3 > #BEGIN OUTFILE4 > #emit mult FULL CONST $eq_no = 1 i = 1 > array_tmp1[kkk] := array_m1[kkk] * array_const_2D0[1]; > #emit mult FULL LINEAR $eq_no = 1 i = 1 > array_tmp2[kkk] := array_tmp1[kkk-1] * array_x[2] + array_tmp1[kkk] * array_x[1]; > #emit mult LINEAR - LINEAR $eq_no = 1 i = 1 > #emit FULL - NOT FULL add $eq_no = 1 > array_tmp4[kkk] := array_tmp3[kkk]; > #emit div FULL FULL $eq_no = 1 > array_tmp5[kkk] := ((array_tmp2[kkk] - ats(kkk,array_tmp4,array_tmp5,2))/array_tmp4[1]); > #emit mult LINEAR - LINEAR $eq_no = 1 i = 1 > #emit FULL - NOT FULL add $eq_no = 1 > array_tmp7[kkk] := array_tmp6[kkk]; > #emit div FULL FULL $eq_no = 1 > array_tmp8[kkk] := ((array_tmp5[kkk] - ats(kkk,array_tmp7,array_tmp8,2))/array_tmp7[1]); > #emit NOT FULL - FULL add $eq_no = 1 > array_tmp9[kkk] := array_tmp8[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_tmp9[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 := 2; > while ((adj2 <= order_d + 1) and (term >= 1)) do # do number 2 > temporary := temporary / glob_h * convfp(adj2); > array_y_higher[adj2,term] := temporary; > adj2 := adj2 + 1; > term := term - 1; > od;# end do number 2 > fi;# end if 2 > fi;# end if 1 > ; > kkk := kkk + 1; > od;# end do number 1 > ; > #BOTTOM ATOMALL > #END OUTFILE4 > #BEGIN OUTFILE5 > > #BOTTOM ATOMALL ??? > # End Function number 9 > end; atomall := proc() local kkk, order_d, adj2, temporary, term; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; array_tmp1[1] := array_m1[1]*array_const_2D0[1]; array_tmp2[1] := array_tmp1[1]*array_x[1]; array_tmp3[1] := array_x[1]*array_x[1]; array_tmp4[1] := array_tmp3[1] + array_const_1D0[1]; array_tmp5[1] := array_tmp2[1]/array_tmp4[1]; array_tmp6[1] := array_x[1]*array_x[1]; array_tmp7[1] := array_tmp6[1] + array_const_1D0[1]; array_tmp8[1] := array_tmp5[1]/array_tmp7[1]; array_tmp9[1] := array_const_0D0[1] + array_tmp8[1]; if not array_y_set_initial[1, 2] then if 1 <= glob_max_terms then temporary := array_tmp9[1]*expt(glob_h, 1)*factorial_3(0, 1); array_y[2] := temporary; array_y_higher[1, 2] := temporary; temporary := temporary*2.0/glob_h; array_y_higher[2, 1] := temporary end if end if; kkk := 2; array_tmp1[2] := array_m1[2]*array_const_2D0[1]; array_tmp2[2] := array_tmp1[2]*array_x[kkk - 1] + array_tmp1[1]*array_x[kkk]; array_tmp3[2] := 2*array_x[1]*array_x[2]; array_tmp4[2] := array_tmp3[2]; array_tmp5[2] := (array_tmp2[2] - ats(2, array_tmp4, array_tmp5, 2))/array_tmp4[1]; array_tmp6[2] := 2*array_x[1]*array_x[2]; array_tmp7[2] := array_tmp6[2]; array_tmp8[2] := (array_tmp5[2] - ats(2, array_tmp7, array_tmp8, 2))/array_tmp7[1]; array_tmp9[2] := array_tmp8[2]; if not array_y_set_initial[1, 3] then if 2 <= glob_max_terms then temporary := array_tmp9[2]*expt(glob_h, 1)*factorial_3(1, 2); array_y[3] := temporary; array_y_higher[1, 3] := temporary; temporary := temporary*2.0/glob_h; array_y_higher[2, 2] := temporary end if end if; kkk := 3; array_tmp1[3] := array_m1[3]*array_const_2D0[1]; array_tmp2[3] := array_tmp1[2]*array_x[kkk - 1] + array_tmp1[1]*array_x[kkk]; array_tmp3[3] := array_x[2]*array_x[2]; array_tmp4[3] := array_tmp3[3]; array_tmp5[3] := (array_tmp2[3] - ats(3, array_tmp4, array_tmp5, 2))/array_tmp4[1]; array_tmp6[3] := array_x[2]*array_x[2]; array_tmp7[3] := array_tmp6[3]; array_tmp8[3] := (array_tmp5[3] - ats(3, array_tmp7, array_tmp8, 2))/array_tmp7[1]; array_tmp9[3] := array_tmp8[3]; if not array_y_set_initial[1, 4] then if 3 <= glob_max_terms then temporary := array_tmp9[3]*expt(glob_h, 1)*factorial_3(2, 3); array_y[4] := temporary; array_y_higher[1, 4] := temporary; temporary := temporary*2.0/glob_h; array_y_higher[2, 3] := temporary end if end if; kkk := 4; array_tmp1[4] := array_m1[4]*array_const_2D0[1]; array_tmp2[4] := array_tmp1[2]*array_x[kkk - 1] + array_tmp1[1]*array_x[kkk]; array_tmp4[4] := array_tmp3[4]; array_tmp5[4] := (array_tmp2[4] - ats(4, array_tmp4, array_tmp5, 2))/array_tmp4[1]; array_tmp7[4] := array_tmp6[4]; array_tmp8[4] := (array_tmp5[4] - ats(4, array_tmp7, array_tmp8, 2))/array_tmp7[1]; array_tmp9[4] := array_tmp8[4]; if not array_y_set_initial[1, 5] then if 4 <= glob_max_terms then temporary := array_tmp9[4]*expt(glob_h, 1)*factorial_3(3, 4); array_y[5] := temporary; array_y_higher[1, 5] := temporary; temporary := temporary*2.0/glob_h; array_y_higher[2, 4] := temporary end if end if; kkk := 5; array_tmp1[5] := array_m1[5]*array_const_2D0[1]; array_tmp2[5] := array_tmp1[2]*array_x[kkk - 1] + array_tmp1[1]*array_x[kkk]; array_tmp4[5] := array_tmp3[5]; array_tmp5[5] := (array_tmp2[5] - ats(5, array_tmp4, array_tmp5, 2))/array_tmp4[1]; array_tmp7[5] := array_tmp6[5]; array_tmp8[5] := (array_tmp5[5] - ats(5, array_tmp7, array_tmp8, 2))/array_tmp7[1]; array_tmp9[5] := array_tmp8[5]; if not array_y_set_initial[1, 6] then if 5 <= glob_max_terms then temporary := array_tmp9[5]*expt(glob_h, 1)*factorial_3(4, 5); array_y[6] := temporary; array_y_higher[1, 6] := temporary; temporary := temporary*2.0/glob_h; array_y_higher[2, 5] := temporary end if end if; kkk := 6; while kkk <= glob_max_terms do array_tmp1[kkk] := array_m1[kkk]*array_const_2D0[1]; array_tmp2[kkk] := array_tmp1[kkk - 1]*array_x[2] + array_tmp1[kkk]*array_x[1]; array_tmp4[kkk] := array_tmp3[kkk]; array_tmp5[kkk] := ( array_tmp2[kkk] - ats(kkk, array_tmp4, array_tmp5, 2))/ array_tmp4[1]; array_tmp7[kkk] := array_tmp6[kkk]; array_tmp8[kkk] := ( array_tmp5[kkk] - ats(kkk, array_tmp7, array_tmp8, 2))/ array_tmp7[1]; array_tmp9[kkk] := array_tmp8[kkk]; order_d := 1; if kkk + order_d + 1 <= glob_max_terms then if not array_y_set_initial[1, kkk + order_d] then temporary := array_tmp9[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 := 2; while adj2 <= order_d + 1 and 1 <= term do temporary := temporary*convfp(adj2)/glob_h; array_y_higher[adj2, term] := temporary; adj2 := adj2 + 1; term := term - 1 end do end if end if; kkk := kkk + 1 end do end proc > #BEGIN ATS LIBRARY BLOCK > omniout_str := proc(iolevel,str) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > printf("%s\n",str); > fi; > # End Function number 1 > end; omniout_str := proc(iolevel, str) global glob_iolevel; if iolevel <= glob_iolevel then printf("%s\n", str) end if end proc > omniout_str_noeol := proc(iolevel,str) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > printf("%s",str); > fi; > # End Function number 1 > end; omniout_str_noeol := proc(iolevel, str) global glob_iolevel; if iolevel <= glob_iolevel then printf("%s", str) end if end proc > omniout_labstr := proc(iolevel,label,str) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > print(label,str); > fi; > # End Function number 1 > end; omniout_labstr := proc(iolevel, label, str) global glob_iolevel; if iolevel <= glob_iolevel then print(label, str) end if end proc > omniout_float := proc(iolevel,prelabel,prelen,value,vallen,postlabel) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > if vallen = 4 then > printf("%-30s = %-42.4g %s \n",prelabel,value, postlabel); > else > printf("%-30s = %-42.32g %s \n",prelabel,value, postlabel); > fi; > fi; > # End Function number 1 > end; omniout_float := proc(iolevel, prelabel, prelen, value, vallen, postlabel) global glob_iolevel; if iolevel <= glob_iolevel then if vallen = 4 then printf("%-30s = %-42.4g %s \n", prelabel, value, postlabel) else printf("%-30s = %-42.32g %s \n", prelabel, value, postlabel) end if end if end proc > omniout_int := proc(iolevel,prelabel,prelen,value,vallen,postlabel) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > if vallen = 5 then > printf("%-30s = %-32d %s\n",prelabel,value, postlabel); > else > printf("%-30s = %-32d %s \n",prelabel,value, postlabel); > fi; > fi; > # End Function number 1 > end; omniout_int := proc(iolevel, prelabel, prelen, value, vallen, postlabel) global glob_iolevel; if iolevel <= glob_iolevel then if vallen = 5 then printf("%-30s = %-32d %s\n", prelabel, value, postlabel) else printf("%-30s = %-32d %s \n", prelabel, value, postlabel) end if end if end proc > omniout_float_arr := proc(iolevel,prelabel,elemnt,prelen,value,vallen,postlabel) > global glob_iolevel; > if (glob_iolevel >= iolevel) then > print(prelabel,"[",elemnt,"]",value, postlabel); > fi; > # End Function number 1 > end; omniout_float_arr := proc( iolevel, prelabel, elemnt, prelen, value, vallen, postlabel) global glob_iolevel; if iolevel <= glob_iolevel then print(prelabel, "[", elemnt, "]", value, postlabel) end if end proc > dump_series := proc(iolevel,dump_label,series_name, > array_series,numb) > global glob_iolevel; > local i; > if (glob_iolevel >= iolevel) then > i := 1; > while (i <= numb) do > print(dump_label,series_name > ,i,array_series[i]); > i := i + 1; > od; > fi; > # End Function number 1 > end; dump_series := proc(iolevel, dump_label, series_name, array_series, numb) local i; global glob_iolevel; if iolevel <= glob_iolevel then i := 1; while i <= numb do print(dump_label, series_name, i, array_series[i]); i := i + 1 end do end if end proc > dump_series_2 := proc(iolevel,dump_label,series_name2, > array_series2,numb,subnum,array_x) > global glob_iolevel; > local i,sub,ts_term; > if (glob_iolevel >= iolevel) then > sub := 1; > while (sub <= subnum) do > i := 1; > while (i <= numb) do > print(dump_label,series_name2,sub,i,array_series2[sub,i]); > od; > sub := sub + 1; > od; > fi; > # End Function number 1 > end; dump_series_2 := proc( iolevel, dump_label, series_name2, array_series2, numb, subnum, array_x) local i, sub, ts_term; global glob_iolevel; if iolevel <= glob_iolevel then sub := 1; while sub <= subnum do i := 1; while i <= numb do print(dump_label, series_name2, sub, i, array_series2[sub, i]) end do; sub := sub + 1 end do end if end proc > cs_info := proc(iolevel,str) > global glob_iolevel,glob_correct_start_flag,glob_h,glob_reached_optimal_h; > if (glob_iolevel >= iolevel) then > print("cs_info " , str , " glob_correct_start_flag = " , glob_correct_start_flag , "glob_h := " , glob_h , "glob_reached_optimal_h := " , glob_reached_optimal_h) > fi; > # End Function number 1 > end; cs_info := proc(iolevel, str) global glob_iolevel, glob_correct_start_flag, glob_h, glob_reached_optimal_h; if iolevel <= glob_iolevel then print("cs_info ", str, " glob_correct_start_flag = ", glob_correct_start_flag, "glob_h := ", glob_h, "glob_reached_optimal_h := ", glob_reached_optimal_h) end if end proc > # Begin Function number 2 > logitem_time := proc(fd,secs_in) > global centuries_in_millinium, days_in_year, hours_in_day, min_in_hour, sec_in_minute, years_in_century; > local cent_int, centuries, days, days_int, hours, hours_int, millinium_int, milliniums, minutes, minutes_int, sec_in_millinium, sec_int, seconds, secs, years, years_int; > secs := (secs_in); > fprintf(fd,""); > if (secs >= 0.0) then # if number 1 > sec_in_millinium := convfloat(sec_in_minute * min_in_hour * hours_in_day * days_in_year * years_in_century * centuries_in_millinium); > milliniums := convfloat(secs / sec_in_millinium); > millinium_int := floor(milliniums); > centuries := (milliniums - millinium_int)*centuries_in_millinium; > cent_int := floor(centuries); > years := (centuries - cent_int) * years_in_century; > years_int := floor(years); > days := (years - years_int) * days_in_year; > days_int := floor(days); > hours := (days - days_int) * hours_in_day; > hours_int := floor(hours); > minutes := (hours - hours_int) * min_in_hour; > minutes_int := floor(minutes); > seconds := (minutes - minutes_int) * sec_in_minute; > sec_int := floor(seconds); > if (millinium_int > 0) then # if number 2 > fprintf(fd,"%d Millinia %d Centuries %d Years %d Days %d Hours %d Minutes %d Seconds",millinium_int,cent_int,years_int,days_int,hours_int,minutes_int,sec_int); > elif (cent_int > 0) then # if number 3 > fprintf(fd,"%d Centuries %d Years %d Days %d Hours %d Minutes %d Seconds",cent_int,years_int,days_int,hours_int,minutes_int,sec_int); > elif (years_int > 0) then # if number 4 > fprintf(fd,"%d Years %d Days %d Hours %d Minutes %d Seconds",years_int,days_int,hours_int,minutes_int,sec_int); > elif (days_int > 0) then # if number 5 > fprintf(fd,"%d Days %d Hours %d Minutes %d Seconds",days_int,hours_int,minutes_int,sec_int); > elif (hours_int > 0) then # if number 6 > fprintf(fd,"%d Hours %d Minutes %d Seconds",hours_int,minutes_int,sec_int); > elif (minutes_int > 0) then # if number 7 > fprintf(fd,"%d Minutes %d Seconds",minutes_int,sec_int); > else > fprintf(fd,"%d Seconds",sec_int); > fi;# end if 7 > else > fprintf(fd,"Unknown"); > fi;# end if 6 > fprintf(fd,""); > # End Function number 2 > end; logitem_time := proc(fd, secs_in) local cent_int, centuries, days, days_int, hours, hours_int, millinium_int, milliniums, minutes, minutes_int, sec_in_millinium, sec_int, seconds, secs, years, years_int; global centuries_in_millinium, days_in_year, hours_in_day, min_in_hour, sec_in_minute, years_in_century; secs := secs_in; fprintf(fd, ""); if 0. <= secs then sec_in_millinium := convfloat(sec_in_minute*min_in_hour* hours_in_day*days_in_year*years_in_century* centuries_in_millinium); milliniums := convfloat(secs/sec_in_millinium); millinium_int := floor(milliniums); centuries := (milliniums - millinium_int)*centuries_in_millinium; cent_int := floor(centuries); years := (centuries - cent_int)*years_in_century; years_int := floor(years); days := (years - years_int)*days_in_year; days_int := floor(days); hours := (days - days_int)*hours_in_day; hours_int := floor(hours); minutes := (hours - hours_int)*min_in_hour; minutes_int := floor(minutes); seconds := (minutes - minutes_int)*sec_in_minute; sec_int := floor(seconds); if 0 < millinium_int then fprintf(fd, "%d Millinia %d Centuries %\ d Years %d Days %d Hours %d Minutes %d Seconds", millinium_int, cent_int, years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < cent_int then fprintf(fd, "%d Centuries %d Years %d Days %d Hours %d Minutes %d Seconds", cent_int, years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < years_int then fprintf(fd, "%d Years %d Days %d Hours %d Minutes %d Seconds", years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < days_int then fprintf(fd, "%d Days %d Hours %d Minutes %d Seconds", days_int, hours_int, minutes_int, sec_int) elif 0 < hours_int then fprintf(fd, "%d Hours %d Minutes %d Seconds", hours_int, minutes_int, sec_int) elif 0 < minutes_int then fprintf(fd, "%d Minutes %d Seconds", minutes_int, sec_int) else fprintf(fd, "%d Seconds", sec_int) end if else fprintf(fd, "Unknown") end if; fprintf(fd, "") end proc > omniout_timestr := proc (secs_in) > global centuries_in_millinium, days_in_year, hours_in_day, min_in_hour, sec_in_minute, years_in_century; > local cent_int, centuries, days, days_int, hours, hours_int, millinium_int, milliniums, minutes, minutes_int, sec_in_millinium, sec_int, seconds, secs, years, years_int; > secs := convfloat(secs_in); > if (secs >= 0.0) then # if number 6 > sec_in_millinium := convfloat(sec_in_minute * min_in_hour * hours_in_day * days_in_year * years_in_century * centuries_in_millinium); > milliniums := convfloat(secs / sec_in_millinium); > millinium_int := floor(milliniums); > centuries := (milliniums - millinium_int)*centuries_in_millinium; > cent_int := floor(centuries); > years := (centuries - cent_int) * years_in_century; > years_int := floor(years); > days := (years - years_int) * days_in_year; > days_int := floor(days); > hours := (days - days_int) * hours_in_day; > hours_int := floor(hours); > minutes := (hours - hours_int) * min_in_hour; > minutes_int := floor(minutes); > seconds := (minutes - minutes_int) * sec_in_minute; > sec_int := floor(seconds); > > if (millinium_int > 0) then # if number 7 > printf(" = %d Millinia %d Centuries %d Years %d Days %d Hours %d Minutes %d Seconds\n",millinium_int,cent_int,years_int,days_int,hours_int,minutes_int,sec_int); > elif (cent_int > 0) then # if number 8 > printf(" = %d Centuries %d Years %d Days %d Hours %d Minutes %d Seconds\n",cent_int,years_int,days_int,hours_int,minutes_int,sec_int); > elif (years_int > 0) then # if number 9 > printf(" = %d Years %d Days %d Hours %d Minutes %d Seconds\n",years_int,days_int,hours_int,minutes_int,sec_int); > elif (days_int > 0) then # if number 10 > printf(" = %d Days %d Hours %d Minutes %d Seconds\n",days_int,hours_int,minutes_int,sec_int); > elif (hours_int > 0) then # if number 11 > printf(" = %d Hours %d Minutes %d Seconds\n",hours_int,minutes_int,sec_int); > elif (minutes_int > 0) then # if number 12 > printf(" = %d Minutes %d Seconds\n",minutes_int,sec_int); > else > printf(" = %d Seconds\n",sec_int); > fi;# end if 12 > else > printf(" Unknown\n"); > fi;# end if 11 > # End Function number 2 > end; omniout_timestr := proc(secs_in) local cent_int, centuries, days, days_int, hours, hours_int, millinium_int, milliniums, minutes, minutes_int, sec_in_millinium, sec_int, seconds, secs, years, years_int; global centuries_in_millinium, days_in_year, hours_in_day, min_in_hour, sec_in_minute, years_in_century; secs := convfloat(secs_in); if 0. <= secs then sec_in_millinium := convfloat(sec_in_minute*min_in_hour* hours_in_day*days_in_year*years_in_century* centuries_in_millinium); milliniums := convfloat(secs/sec_in_millinium); millinium_int := floor(milliniums); centuries := (milliniums - millinium_int)*centuries_in_millinium; cent_int := floor(centuries); years := (centuries - cent_int)*years_in_century; years_int := floor(years); days := (years - years_int)*days_in_year; days_int := floor(days); hours := (days - days_int)*hours_in_day; hours_int := floor(hours); minutes := (hours - hours_int)*min_in_hour; minutes_int := floor(minutes); seconds := (minutes - minutes_int)*sec_in_minute; sec_int := floor(seconds); if 0 < millinium_int then printf(" = %d Millinia %d Centuries %d\ Years %d Days %d Hours %d Minutes %d Seconds\n", millinium_int, cent_int, years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < cent_int then printf(" = %d Centuries %d Years %d Days \ %d Hours %d Minutes %d Seconds\n", cent_int, years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < years_int then printf( " = %d Years %d Days %d Hours %d Minutes %d Seconds\n", years_int, days_int, hours_int, minutes_int, sec_int) elif 0 < days_int then printf( " = %d Days %d Hours %d Minutes %d Seconds\n", days_int, hours_int, minutes_int, sec_int) elif 0 < hours_int then printf( " = %d Hours %d Minutes %d Seconds\n", hours_int, minutes_int, sec_int) elif 0 < minutes_int then printf(" = %d Minutes %d Seconds\n", minutes_int, sec_int) else printf(" = %d Seconds\n", sec_int) end if else printf(" Unknown\n") end if end proc > # Begin Function number 3 > ats := proc( > mmm_ats,array_a,array_b,jjj_ats) > local iii_ats, lll_ats,ma_ats, ret_ats; > > > > > > ret_ats := 0.0; > if (jjj_ats <= mmm_ats) then # if number 11 > 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 + array_a[iii_ats]*array_b[lll_ats]; > iii_ats := iii_ats + 1; > od;# end do number 1 > fi;# end if 11 > ; > ret_ats; > > # End Function number 3 > end; ats := proc(mmm_ats, array_a, array_b, jjj_ats) local iii_ats, lll_ats, ma_ats, ret_ats; ret_ats := 0.; if jjj_ats <= mmm_ats then ma_ats := mmm_ats + 1; iii_ats := jjj_ats; while iii_ats <= mmm_ats do lll_ats := ma_ats - iii_ats; ret_ats := ret_ats + array_a[iii_ats]*array_b[lll_ats]; iii_ats := iii_ats + 1 end do end if; ret_ats end proc > # Begin Function number 4 > att := proc( > mmm_att,array_aa,array_bb,jjj_att) > global glob_max_terms; > local al_att, iii_att,lll_att, ma_att, ret_att; > > > > > > ret_att := 0.0; > if (jjj_att <= mmm_att) then # if number 11 > 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 12 > ret_att := ret_att + array_aa[iii_att]*array_bb[lll_att]* convfp(al_att); > fi;# end if 12 > ; > iii_att := iii_att + 1; > od;# end do number 1 > ; > ret_att := ret_att / convfp(mmm_att) ; > fi;# end if 11 > ; > ret_att; > > # End Function number 4 > end; att := proc(mmm_att, array_aa, array_bb, jjj_att) local al_att, iii_att, lll_att, ma_att, ret_att; global glob_max_terms; ret_att := 0.; if jjj_att <= mmm_att then ma_att := mmm_att + 2; iii_att := jjj_att; while iii_att <= mmm_att do lll_att := ma_att - iii_att; al_att := lll_att - 1; if lll_att <= glob_max_terms then ret_att := ret_att + array_aa[iii_att]*array_bb[lll_att]*convfp(al_att) end if; iii_att := iii_att + 1 end do; ret_att := ret_att/convfp(mmm_att) end if; ret_att end proc > # Begin Function number 5 > display_pole := proc() > global ALWAYS,glob_display_flag, glob_large_float, array_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 11 > omniout_float(ALWAYS,"Radius of convergence ",4, array_pole[1],4," "); > omniout_float(ALWAYS,"Order of pole ",4, array_pole[2],4," "); > fi;# end if 11 > # End Function number 5 > end; display_pole := proc() global ALWAYS, glob_display_flag, glob_large_float, array_pole; if array_pole[1] <> glob_large_float and 0. < array_pole[1] and array_pole[2] <> glob_large_float and 0. < array_pole[2] and glob_display_flag then omniout_float(ALWAYS, "Radius of convergence ", 4, array_pole[1], 4, " "); omniout_float(ALWAYS, "Order of pole ", 4, array_pole[2], 4, " ") end if end proc > # Begin Function number 6 > logditto := proc(file) > fprintf(file,""); > fprintf(file,"ditto"); > fprintf(file,""); > # End Function number 6 > end; logditto := proc(file) fprintf(file, ""); fprintf(file, "ditto"); fprintf(file, "") end proc > # Begin Function number 7 > logitem_integer := proc(file,n) > fprintf(file,""); > fprintf(file,"%d",n); > fprintf(file,""); > # End Function number 7 > end; logitem_integer := proc(file, n) fprintf(file, ""); fprintf(file, "%d", n); fprintf(file, "") end proc > # Begin Function number 8 > logitem_str := proc(file,str) > fprintf(file,""); > fprintf(file,str); > fprintf(file,""); > # End Function number 8 > end; logitem_str := proc(file, str) fprintf(file, ""); fprintf(file, str); fprintf(file, "") end proc > # Begin Function number 9 > logitem_good_digits := proc(file,rel_error) > global glob_small_float; > > local good_digits; > > > fprintf(file,""); > if (rel_error <> -1.0) then # if number 11 > if (rel_error <> 0.0) then # if number 12 > good_digits := -trunc(log10(rel_error/100.0)); > fprintf(file,"%d",good_digits); > else > good_digits := Digits; > fprintf(file,"%d",good_digits); > fi;# end if 12 > ; > else > fprintf(file,"Unknown"); > fi;# end if 11 > ; > fprintf(file,""); > > # End Function number 9 > end; logitem_good_digits := proc(file, rel_error) local good_digits; global glob_small_float; fprintf(file, ""); if rel_error <> -1.0 then if rel_error <> 0. then good_digits := -trunc(log10(rel_error/100.0)); fprintf(file, "%d", good_digits) else good_digits := Digits; fprintf(file, "%d", good_digits) end if else fprintf(file, "Unknown") end if; fprintf(file, "") end proc > # Begin Function number 10 > log_revs := proc(file,revs) > fprintf(file,revs); > # End Function number 10 > end; log_revs := proc(file, revs) fprintf(file, revs) end proc > # Begin Function number 11 > logitem_float := proc(file,x) > fprintf(file,""); > fprintf(file,"%g",x); > fprintf(file,""); > # End Function number 11 > end; logitem_float := proc(file, x) fprintf(file, ""); fprintf(file, "%g", x); fprintf(file, "") end proc > # Begin Function number 12 > logitem_pole := proc(file,pole) > fprintf(file,""); > if (pole = 0) then # if number 11 > fprintf(file,"NA"); > elif (pole = 1) then # if number 12 > fprintf(file,"Real"); > elif (pole = 2) then # if number 13 > fprintf(file,"Complex"); > else > fprintf(file,"No Pole"); > fi;# end if 13 > fprintf(file,""); > # End Function number 12 > end; logitem_pole := proc(file, pole) fprintf(file, ""); if pole = 0 then fprintf(file, "NA") elif pole = 1 then fprintf(file, "Real") elif pole = 2 then fprintf(file, "Complex") else fprintf(file, "No Pole") end if; fprintf(file, "") end proc > # Begin Function number 13 > logstart := proc(file) > fprintf(file,""); > # End Function number 13 > end; logstart := proc(file) fprintf(file, "") end proc > # Begin Function number 14 > logend := proc(file) > fprintf(file,"\n"); > # End Function number 14 > end; logend := proc(file) fprintf(file, "\n") end proc > # Begin Function number 15 > not_reached_end := proc(x,x_end) > global neg_h; > local ret; > > > > if ((glob_neg_h and (x > x_end)) or (( not glob_neg_h) and (x < x_end))) then # if number 13 > ret := true; > else > ret := false; > fi;# end if 13 > ; > > > ret; > > # End Function number 15 > end; not_reached_end := proc(x, x_end) local ret; global neg_h; if glob_neg_h and x_end < x or not glob_neg_h and x < x_end then ret := true else ret := false end if; ret end proc > # Begin Function number 16 > chk_data := proc() > global glob_max_iter,ALWAYS, glob_max_terms; > local errflag; > > > > errflag := false; > > if ((glob_max_terms < 15) or (glob_max_terms > 512)) then # if number 13 > omniout_str(ALWAYS,"Illegal max_terms = -- Using 30"); > glob_max_terms := 30; > fi;# end if 13 > ; > if (glob_max_iter < 2) then # if number 13 > omniout_str(ALWAYS,"Illegal max_iter"); > errflag := true; > fi;# end if 13 > ; > if (errflag) then # if number 13 > > quit; > fi;# end if 13 > > # End Function number 16 > end; chk_data := proc() local errflag; global glob_max_iter, ALWAYS, glob_max_terms; errflag := false; if glob_max_terms < 15 or 512 < glob_max_terms then omniout_str(ALWAYS, "Illegal max_terms = -- Using 30"); glob_max_terms := 30 end if; if glob_max_iter < 2 then omniout_str(ALWAYS, "Illegal max_iter"); errflag := true end if; if errflag then quit end if end proc > # Begin Function number 17 > comp_expect_sec := proc(t_end2,t_start2,t2,clock_sec2) > global glob_small_float; > local ms2, rrr, sec_left, sub1, sub2; > > > > ; > ms2 := clock_sec2; > sub1 := (t_end2-t_start2); > sub2 := (t2-t_start2); > if (sub1 = 0.0) then # if number 13 > sec_left := 0.0; > else > if (sub2 > 0.0) then # if number 14 > rrr := (sub1/sub2); > sec_left := rrr * ms2 - ms2; > else > sec_left := 0.0; > fi;# end if 14 > fi;# end if 13 > ; > sec_left; > > # End Function number 17 > end; comp_expect_sec := proc(t_end2, t_start2, t2, clock_sec2) local ms2, rrr, sec_left, sub1, sub2; global glob_small_float; ms2 := clock_sec2; sub1 := t_end2 - t_start2; sub2 := t2 - t_start2; if sub1 = 0. then sec_left := 0. else if 0. < sub2 then rrr := sub1/sub2; sec_left := rrr*ms2 - ms2 else sec_left := 0. end if end if; sec_left end proc > # Begin Function number 18 > comp_percent := proc(t_end2,t_start2, t2) > global glob_small_float; > local rrr, sub1, sub2; > > > > sub1 := (t_end2-t_start2); > sub2 := (t2-t_start2); > if (sub2 > glob_small_float) then # if number 13 > rrr := (100.0*sub2)/sub1; > else > rrr := 0.0; > fi;# end if 13 > ; > rrr; > > # End Function number 18 > end; comp_percent := proc(t_end2, t_start2, t2) local rrr, sub1, sub2; global glob_small_float; sub1 := t_end2 - t_start2; sub2 := t2 - t_start2; if glob_small_float < sub2 then rrr := 100.0*sub2/sub1 else rrr := 0. end if; rrr end proc > # Begin Function number 19 > factorial_2 := proc(nnn) > local ret; > > > > ret := nnn!; > > # End Function number 19 > end; factorial_2 := proc(nnn) local ret; ret := nnn! end proc > # Begin Function number 20 > factorial_1 := proc(nnn) > global glob_max_terms,array_fact_1; > local ret; > > > > if (nnn <= glob_max_terms) then # if number 13 > if (array_fact_1[nnn] = 0) then # if number 14 > ret := factorial_2(nnn); > array_fact_1[nnn] := ret; > else > ret := array_fact_1[nnn]; > fi;# end if 14 > ; > else > ret := factorial_2(nnn); > fi;# end if 13 > ; > ret; > > # End Function number 20 > end; factorial_1 := proc(nnn) local ret; global glob_max_terms, array_fact_1; if nnn <= glob_max_terms then if array_fact_1[nnn] = 0 then ret := factorial_2(nnn); array_fact_1[nnn] := ret else ret := array_fact_1[nnn] end if else ret := factorial_2(nnn) end if; ret end proc > # Begin Function number 21 > factorial_3 := proc(mmm,nnn) > global glob_max_terms,array_fact_2; > local ret; > > > > if ((nnn <= glob_max_terms) and (mmm <= glob_max_terms)) then # if number 13 > if (array_fact_2[mmm,nnn] = 0) then # if number 14 > ret := factorial_1(mmm)/factorial_1(nnn); > array_fact_2[mmm,nnn] := ret; > else > ret := array_fact_2[mmm,nnn]; > fi;# end if 14 > ; > else > ret := factorial_2(mmm)/factorial_2(nnn); > fi;# end if 13 > ; > ret; > > # End Function number 21 > end; factorial_3 := proc(mmm, nnn) local ret; global glob_max_terms, array_fact_2; if nnn <= glob_max_terms and mmm <= glob_max_terms then if array_fact_2[mmm, nnn] = 0 then ret := factorial_1(mmm)/factorial_1(nnn); array_fact_2[mmm, nnn] := ret else ret := array_fact_2[mmm, nnn] end if else ret := factorial_2(mmm)/factorial_2(nnn) end if; ret end proc > # Begin Function number 22 > convfp := proc(mmm) > (mmm); > > # End Function number 22 > end; convfp := proc(mmm) mmm end proc > # Begin Function number 23 > convfloat := proc(mmm) > (mmm); > > # End Function number 23 > end; convfloat := proc(mmm) mmm end proc > elapsed_time_seconds := proc() > time(); > end; elapsed_time_seconds := proc() time() end proc > > > > > > > > > > > > > > > > > > > > omniabs := proc(x) > abs(x); > end; omniabs := proc(x) abs(x) end proc > expt := proc(x,y) > (x^y); > end; expt := proc(x, y) x^y end proc > #END ATS LIBRARY BLOCK > #BEGIN USER DEF BLOCK > #BEGIN USER DEF BLOCK > exact_soln_y := proc(x) > return(1.0 / (x * x + 1.0)); > end; exact_soln_y := proc(x) return 1.0/(x*x + 1.0) end proc > > > #END USER DEF BLOCK > #END USER DEF BLOCK > #END OUTFILE5 > # Begin Function number 2 > main := proc() > #BEGIN OUTFIEMAIN > local d1,d2,d3,d4,est_err_2,niii,done_once, > term,ord,order_diff,term_no,html_log_file,iiif,jjjf, > rows,r_order,sub_iter,calc_term,iii,temp_sum,current_iter, > x_start,x_end > ,it, log10norm, max_terms, opt_iter, tmp,subiter; > global > ALWAYS, > DEBUGL, > DEBUGMASSIVE, > INFO, > glob_max_terms, > glob_iolevel, > #Top Generate Globals Decl > glob_log10normmin, > glob_normmax, > glob_relerr, > glob_abserr, > glob_optimal_done, > glob_warned, > glob_max_rel_trunc_err, > glob_log10_relerr, > glob_neg_h, > glob_hmax, > glob_not_yet_start_msg, > glob_max_minutes, > glob_curr_iter_when_opt, > glob_optimal_clock_start_sec, > glob_next_display, > glob_look_poles, > glob_last_good_h, > glob_disp_incr, > glob_current_iter, > glob_smallish_float, > glob_max_trunc_err, > glob_max_iter, > glob_max_hours, > glob_not_yet_finished, > centuries_in_millinium, > djd_debug, > glob_max_sec, > glob_warned2, > glob_no_eqs, > glob_log10_abserr, > sec_in_minute, > glob_max_opt_iter, > glob_optimal_expect_sec, > glob_start, > glob_dump_analytic, > days_in_year, > MAX_UNCHANGED, > glob_optimal_start, > glob_large_float, > glob_clock_start_sec, > djd_debug2, > glob_good_digits, > glob_html_log, > glob_percent_done, > glob_log10abserr, > glob_iter, > hours_in_day, > glob_subiter_method, > glob_orig_start_sec, > glob_small_float, > glob_reached_optimal_h, > years_in_century, > glob_display_flag, > glob_log10relerr, > glob_unchanged_h_cnt, > glob_display_interval, > glob_hmin_init, > glob_hmin, > glob_h, > glob_initial_pass, > glob_clock_sec, > glob_almost_1, > min_in_hour, > glob_dump, > #Bottom Generate Globals Decl > #BEGIN CONST > array_const_1D0, > array_const_1, > array_const_2D0, > array_const_0D0, > #END CONST > array_y, > array_x, > array_norms, > array_y_init, > array_last_rel_error, > array_tmp0, > array_tmp1, > array_tmp2, > array_tmp3, > array_tmp4, > array_tmp5, > array_tmp6, > array_tmp7, > array_tmp8, > array_tmp9, > array_type_pole, > array_m1, > array_fact_1, > array_pole, > array_1st_rel_error, > array_y_higher, > array_real_pole, > array_complex_pole, > array_y_higher_work2, > array_fact_2, > array_poles, > array_y_set_initial, > array_y_higher_work, > glob_last; > glob_last; > ALWAYS := 1; > INFO := 2; > DEBUGL := 3; > DEBUGMASSIVE := 4; > glob_iolevel := INFO; > ALWAYS := 1; > DEBUGL := 3; > DEBUGMASSIVE := 4; > INFO := 2; > glob_max_terms := 30; > glob_iolevel := 5; > glob_log10normmin := 0.1; > glob_normmax := 0.0; > glob_relerr := 0.1e-10; > glob_abserr := 0.1e-10; > glob_optimal_done := false; > glob_warned := false; > glob_max_rel_trunc_err := 0.1e-10; > glob_log10_relerr := 0.1e-10; > glob_neg_h := false; > glob_hmax := 1.0; > glob_not_yet_start_msg := true; > glob_max_minutes := 0.0; > glob_curr_iter_when_opt := 0; > glob_optimal_clock_start_sec := 0.0; > glob_next_display := 0.0; > glob_look_poles := false; > glob_last_good_h := 0.1; > glob_disp_incr := 0.1; > glob_current_iter := 0; > glob_smallish_float := 0.1e-100; > glob_max_trunc_err := 0.1e-10; > glob_max_iter := 1000; > glob_max_hours := 0.0; > glob_not_yet_finished := true; > centuries_in_millinium := 10; > djd_debug := true; > glob_max_sec := 10000.0; > glob_warned2 := false; > glob_no_eqs := 0; > glob_log10_abserr := 0.1e-10; > sec_in_minute := 60; > glob_max_opt_iter := 10; > glob_optimal_expect_sec := 0.1; > glob_start := 0; > glob_dump_analytic := false; > days_in_year := 365; > MAX_UNCHANGED := 10; > glob_optimal_start := 0.0; > glob_large_float := 9.0e100; > glob_clock_start_sec := 0.0; > djd_debug2 := true; > glob_good_digits := 0; > glob_html_log := true; > glob_percent_done := 0.0; > glob_log10abserr := 0.0; > glob_iter := 0; > hours_in_day := 24; > glob_subiter_method := 3; > glob_orig_start_sec := 0.0; > glob_small_float := 0.1e-50; > glob_reached_optimal_h := false; > years_in_century := 100; > glob_display_flag := true; > glob_log10relerr := 0.0; > glob_unchanged_h_cnt := 0; > glob_display_interval := 0.0; > glob_hmin_init := 0.001; > glob_hmin := 0.00000000001; > glob_h := 0.1; > glob_initial_pass := true; > glob_clock_sec := 0.0; > glob_almost_1 := 0.9990; > min_in_hour := 60; > glob_dump := false; > #Write Set Defaults > glob_orig_start_sec := elapsed_time_seconds(); > MAX_UNCHANGED := 10; > 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/sing4postode.ode#################"); > omniout_str(ALWAYS,"diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0) /( x * x + 1.0);"); > omniout_str(ALWAYS,"!"); > omniout_str(ALWAYS,"#BEGIN FIRST INPUT BLOCK"); > omniout_str(ALWAYS,"Digits := 32;"); > 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 := -2.0;"); > omniout_str(ALWAYS,"x_end := 1.0;"); > omniout_str(ALWAYS,"array_y_init[0 + 1] := exact_soln_y(x_start);"); > omniout_str(ALWAYS,"glob_h := 0.1;"); > omniout_str(ALWAYS,"glob_look_poles := true;"); > omniout_str(ALWAYS,"glob_max_iter := 50;"); > omniout_str(ALWAYS,"#END SECOND INPUT BLOCK"); > omniout_str(ALWAYS,"#BEGIN OVERRIDE BLOCK"); > omniout_str(ALWAYS,"glob_h := 0.005 ;"); > omniout_str(ALWAYS,"glob_display_interval := 0.1;"); > omniout_str(ALWAYS,"glob_look_poles := true;"); > omniout_str(ALWAYS,"glob_max_iter := 10000;"); > omniout_str(ALWAYS,"glob_max_minutes := 10;"); > omniout_str(ALWAYS,"#END OVERRIDE BLOCK"); > omniout_str(ALWAYS,"!"); > omniout_str(ALWAYS,"#BEGIN USER DEF BLOCK"); > omniout_str(ALWAYS,"exact_soln_y := proc(x)"); > omniout_str(ALWAYS,"return(1.0 / (x * x + 1.0));"); > omniout_str(ALWAYS,"end;"); > omniout_str(ALWAYS,""); > omniout_str(ALWAYS,""); > omniout_str(ALWAYS,"#END USER DEF BLOCK"); > omniout_str(ALWAYS,"#######END OF ECHO OF PROBLEM#################"); > glob_unchanged_h_cnt := 0; > glob_warned := false; > glob_warned2 := false; > glob_small_float := 1.0e-200; > glob_smallish_float := 1.0e-64; > glob_large_float := 1.0e100; > glob_almost_1 := 0.99; > glob_log10_abserr := -8.0; > glob_log10_relerr := -8.0; > glob_hmax := 0.01; > #BEGIN FIRST INPUT BLOCK > #BEGIN FIRST INPUT BLOCK > Digits := 32; > 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:= Array(0..(max_terms + 1),[]); > array_x:= Array(0..(max_terms + 1),[]); > array_norms:= Array(0..(max_terms + 1),[]); > array_y_init:= Array(0..(max_terms + 1),[]); > array_last_rel_error:= Array(0..(max_terms + 1),[]); > array_tmp0:= Array(0..(max_terms + 1),[]); > array_tmp1:= Array(0..(max_terms + 1),[]); > array_tmp2:= Array(0..(max_terms + 1),[]); > array_tmp3:= Array(0..(max_terms + 1),[]); > array_tmp4:= Array(0..(max_terms + 1),[]); > array_tmp5:= Array(0..(max_terms + 1),[]); > array_tmp6:= Array(0..(max_terms + 1),[]); > array_tmp7:= Array(0..(max_terms + 1),[]); > array_tmp8:= Array(0..(max_terms + 1),[]); > array_tmp9:= Array(0..(max_terms + 1),[]); > array_type_pole:= Array(0..(max_terms + 1),[]); > array_m1:= Array(0..(max_terms + 1),[]); > array_fact_1:= Array(0..(max_terms + 1),[]); > array_pole:= Array(0..(max_terms + 1),[]); > array_1st_rel_error:= Array(0..(max_terms + 1),[]); > array_y_higher := Array(0..(2+ 1) ,(0..max_terms+ 1),[]); > array_real_pole := Array(0..(1+ 1) ,(0..3+ 1),[]); > array_complex_pole := Array(0..(1+ 1) ,(0..3+ 1),[]); > array_y_higher_work2 := Array(0..(2+ 1) ,(0..max_terms+ 1),[]); > array_fact_2 := Array(0..(max_terms+ 1) ,(0..max_terms+ 1),[]); > array_poles := Array(0..(1+ 1) ,(0..3+ 1),[]); > array_y_set_initial := Array(0..(2+ 1) ,(0..max_terms+ 1),[]); > array_y_higher_work := Array(0..(2+ 1) ,(0..max_terms+ 1),[]); > term := 1; > while (term <= max_terms) do # do number 2 > array_y[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_x[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_norms[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_y_init[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_last_rel_error[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp0[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp2[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp3[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp4[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp5[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp6[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp7[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp8[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_tmp9[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_type_pole[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_m1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_fact_1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_pole[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > term := 1; > while (term <= max_terms) do # do number 2 > array_1st_rel_error[term] := 0.0; > term := term + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# 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; > od;# end do number 3 > ; > ord := ord + 1; > od;# end do number 2 > ; > #BEGIN ARRAYS DEFINED AND INITIALIZATED > array_x := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_x[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_y := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_y[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp9 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp9[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp8 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp8[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp7 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp7[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp6 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp6[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp5 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp5[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp4 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp4[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp3 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp3[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp2 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp2[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp1 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_tmp0 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_tmp0[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_m1 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_m1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_const_1D0 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_const_1D0[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_const_1D0[1] := 1.0; > array_const_1 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_const_1[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_const_1[1] := 1; > array_const_2D0 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_const_2D0[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_const_2D0[1] := 2.0; > array_const_0D0 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms + 1) do # do number 2 > array_const_0D0[term] := 0.0; > term := term + 1; > od;# end do number 2 > ; > array_const_0D0[1] := 0.0; > array_m1 := Array(1..(max_terms+1 + 1),[]); > term := 1; > while (term <= max_terms) do # do number 2 > array_m1[term] := 0.0; > term := term + 1; > od;# 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; > od;# end do number 3 > ; > iiif := iiif + 1; > od;# end do number 2 > ; > #Done Initing Factorial Tables > #TOP SECOND INPUT BLOCK > #BEGIN SECOND INPUT BLOCK > #END FIRST INPUT BLOCK > #BEGIN SECOND INPUT BLOCK > x_start := -2.0; > x_end := 1.0; > array_y_init[0 + 1] := exact_soln_y(x_start); > glob_h := 0.1; > glob_look_poles := true; > glob_max_iter := 50; > #END SECOND INPUT BLOCK > #BEGIN OVERRIDE BLOCK > glob_h := 0.005 ; > glob_display_interval := 0.1; > glob_look_poles := true; > glob_max_iter := 10000; > glob_max_minutes := 10; > #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); > fi;# 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; > if (glob_html_log) then # if number 2 > html_log_file := fopen("html/entry.html",WRITE,TEXT); > fi;# end if 2 > ; > #BEGIN SOLUTION CODE > 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; > od;# 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; > od;# end do number 3 > ; > r_order := r_order + 1; > od;# end do number 2 > ; > current_iter := 1; > glob_clock_start_sec := elapsed_time_seconds(); > if (omniabs(array_y_higher[1,1]) > glob_small_float) then # if number 2 > tmp := omniabs(array_y_higher[1,1]); > log10norm := (log10(tmp)); > if (log10norm < glob_log10normmin) then # if number 3 > glob_log10normmin := log10norm; > fi;# end if 3 > fi;# end if 2 > ; > 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 not_reached_end(array_x[1] , 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 2 > omniout_str(INFO," "); > omniout_str(INFO,"TOP MAIN SOLVE Loop"); > fi;# end if 2 > ; > glob_iter := glob_iter + 1; > glob_clock_sec := elapsed_time_seconds(); > glob_current_iter := glob_current_iter + 1; > atomall(); > if (glob_look_poles) then # if number 2 > #left paren 0004C > check_for_pole(); > fi;# end if 2 > ;#was right paren 0004C > if (reached_interval()) then # if number 2 > glob_next_display := glob_next_display + glob_display_interval; > fi;# end if 2 > ; > array_x[1] := array_x[1] + glob_h; > array_x[2] := glob_h; > #Jump Series array_y > order_diff := 1; > #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_subseriesarray_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; > od;# end do number 3 > ; > #AFTER ADJUST SUBSERIES EQ =1 > #BEFORE SUM SUBSERIES EQ =1 > temp_sum := 0.0; > ord := 2; > calc_term := 1; > #sum_subseriesarray_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; > od;# 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_subseriesarray_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; > od;# end do number 3 > ; > #AFTER ADJUST SUBSERIES EQ =1 > #BEFORE SUM SUBSERIES EQ =1 > temp_sum := 0.0; > ord := 1; > calc_term := 2; > #sum_subseriesarray_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; > od;# 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_subseriesarray_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; > od;# end do number 3 > ; > #AFTER ADJUST SUBSERIES EQ =1 > #BEFORE SUM SUBSERIES EQ =1 > temp_sum := 0.0; > ord := 1; > calc_term := 1; > #sum_subseriesarray_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; > od;# 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; > od;# end do number 4 > ; > term_no := term_no - 1; > od;# end do number 3 > ; > #END PART 2 HEVE MOVED TERMS to REGULAR Array > display_alot(current_iter) > ; > od;# end do number 2 > ;#right paren 0001C > omniout_str(ALWAYS,"Finished!"); > if (glob_iter >= glob_max_iter) then # if number 2 > omniout_str(ALWAYS,"Maximum Iterations Reached before Solution Completed!"); > fi;# end if 2 > ; > if (elapsed_time_seconds() - convfloat(glob_orig_start_sec) >= convfloat(glob_max_sec )) then # if number 2 > omniout_str(ALWAYS,"Maximum Time Reached before Solution Completed!"); > fi;# end if 2 > ; > glob_clock_sec := elapsed_time_seconds(); > omniout_str(INFO,"diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0) /( x * x + 1.0);"); > omniout_int(INFO,"Iterations ",32,glob_iter,4," ") > ; > prog_report(x_start,x_end); > if (glob_html_log) then # if number 2 > logstart(html_log_file); > logitem_str(html_log_file,"2012-09-21T03:44:08-05:00") > ; > logitem_str(html_log_file,"Maple") > ; > logitem_str(html_log_file,"sing4") > ; > logitem_str(html_log_file,"diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0) /( x * x + 1.0);") > ; > logitem_float(html_log_file,x_start) > ; > logitem_float(html_log_file,x_end) > ; > logitem_float(html_log_file,array_x[1]) > ; > logitem_float(html_log_file,glob_h) > ; > logitem_integer(html_log_file,Digits) > ; > ; > 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 3 > 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; > fi;# end if 3 > ; > logitem_time(html_log_file,convfloat(glob_clock_sec)) > ; > if (glob_percent_done < 100.0) then # if number 3 > logitem_time(html_log_file,convfloat(glob_optimal_expect_sec)) > ; > 0; > else > logitem_str(html_log_file,"Done") > ; > 0; > fi;# end if 3 > ; > log_revs(html_log_file," 130 ") > ; > logitem_str(html_log_file,"sing4 diffeq.mxt") > ; > logitem_str(html_log_file,"sing4 maple results") > ; > logitem_str(html_log_file,"c c++ Maple and Maxima") > ; > logend(html_log_file) > ; > ; > fi;# end if 2 > ; > if (glob_html_log) then # if number 2 > fclose(html_log_file); > fi;# end if 2 > ; > ;; > #END OUTFILEMAIN > > # End Function number 9 > end; main := proc() local d1, d2, d3, d4, est_err_2, niii, done_once, term, ord, order_diff, term_no, html_log_file, iiif, jjjf, rows, r_order, sub_iter, calc_term, iii, temp_sum, current_iter, x_start, x_end, it, log10norm, max_terms, opt_iter, tmp, subiter; global ALWAYS, DEBUGL, DEBUGMASSIVE, INFO, glob_max_terms, glob_iolevel, glob_log10normmin, glob_normmax, glob_relerr, glob_abserr, glob_optimal_done, glob_warned, glob_max_rel_trunc_err, glob_log10_relerr, glob_neg_h, glob_hmax, glob_not_yet_start_msg, glob_max_minutes, glob_curr_iter_when_opt, glob_optimal_clock_start_sec, glob_next_display, glob_look_poles, glob_last_good_h, glob_disp_incr, glob_current_iter, glob_smallish_float, glob_max_trunc_err, glob_max_iter, glob_max_hours, glob_not_yet_finished, centuries_in_millinium, djd_debug, glob_max_sec, glob_warned2, glob_no_eqs, glob_log10_abserr, sec_in_minute, glob_max_opt_iter, glob_optimal_expect_sec, glob_start, glob_dump_analytic, days_in_year, MAX_UNCHANGED, glob_optimal_start, glob_large_float, glob_clock_start_sec, djd_debug2, glob_good_digits, glob_html_log, glob_percent_done, glob_log10abserr, glob_iter, hours_in_day, glob_subiter_method, glob_orig_start_sec, glob_small_float, glob_reached_optimal_h, years_in_century, glob_display_flag, glob_log10relerr, glob_unchanged_h_cnt, glob_display_interval, glob_hmin_init, glob_hmin, glob_h, glob_initial_pass, glob_clock_sec, glob_almost_1, min_in_hour, glob_dump, array_const_1D0, array_const_1, array_const_2D0, array_const_0D0, array_y, array_x, array_norms, array_y_init, array_last_rel_error, array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5, array_tmp6, array_tmp7, array_tmp8, array_tmp9, array_type_pole, array_m1, array_fact_1, array_pole, array_1st_rel_error, array_y_higher, array_real_pole, array_complex_pole, array_y_higher_work2, array_fact_2, array_poles, array_y_set_initial, array_y_higher_work, glob_last; glob_last; ALWAYS := 1; INFO := 2; DEBUGL := 3; DEBUGMASSIVE := 4; glob_iolevel := INFO; ALWAYS := 1; DEBUGL := 3; DEBUGMASSIVE := 4; INFO := 2; glob_max_terms := 30; glob_iolevel := 5; glob_log10normmin := 0.1; glob_normmax := 0.; glob_relerr := 0.1*10^(-10); glob_abserr := 0.1*10^(-10); glob_optimal_done := false; glob_warned := false; glob_max_rel_trunc_err := 0.1*10^(-10); glob_log10_relerr := 0.1*10^(-10); glob_neg_h := false; glob_hmax := 1.0; glob_not_yet_start_msg := true; glob_max_minutes := 0.; glob_curr_iter_when_opt := 0; glob_optimal_clock_start_sec := 0.; glob_next_display := 0.; glob_look_poles := false; glob_last_good_h := 0.1; glob_disp_incr := 0.1; glob_current_iter := 0; glob_smallish_float := 0.1*10^(-100); glob_max_trunc_err := 0.1*10^(-10); glob_max_iter := 1000; glob_max_hours := 0.; glob_not_yet_finished := true; centuries_in_millinium := 10; djd_debug := true; glob_max_sec := 10000.0; glob_warned2 := false; glob_no_eqs := 0; glob_log10_abserr := 0.1*10^(-10); sec_in_minute := 60; glob_max_opt_iter := 10; glob_optimal_expect_sec := 0.1; glob_start := 0; glob_dump_analytic := false; days_in_year := 365; MAX_UNCHANGED := 10; glob_optimal_start := 0.; glob_large_float := 0.90*10^101; glob_clock_start_sec := 0.; djd_debug2 := true; glob_good_digits := 0; glob_html_log := true; glob_percent_done := 0.; glob_log10abserr := 0.; glob_iter := 0; hours_in_day := 24; glob_subiter_method := 3; glob_orig_start_sec := 0.; glob_small_float := 0.1*10^(-50); glob_reached_optimal_h := false; years_in_century := 100; glob_display_flag := true; glob_log10relerr := 0.; glob_unchanged_h_cnt := 0; glob_display_interval := 0.; glob_hmin_init := 0.001; glob_hmin := 0.1*10^(-10); glob_h := 0.1; glob_initial_pass := true; glob_clock_sec := 0.; glob_almost_1 := 0.9990; min_in_hour := 60; glob_dump := false; glob_orig_start_sec := elapsed_time_seconds(); MAX_UNCHANGED := 10; 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.; glob_max_minutes := 15.0; omniout_str(ALWAYS, "##############ECHO OF PROBLEM#################"); omniout_str(ALWAYS, "##############temp/sing4postode.ode#################"); omniout_str(ALWAYS, "diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.\ 0) /( x * x + 1.0);"); omniout_str(ALWAYS, "!"); omniout_str(ALWAYS, "#BEGIN FIRST INPUT BLOCK"); omniout_str(ALWAYS, "Digits := 32;"); 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 := -2.0;"); omniout_str(ALWAYS, "x_end := 1.0;"); omniout_str(ALWAYS, "array_y_init[0 + 1] := exact_soln_y(x_start);"); omniout_str(ALWAYS, "glob_h := 0.1;"); omniout_str(ALWAYS, "glob_look_poles := true;"); omniout_str(ALWAYS, "glob_max_iter := 50;"); omniout_str(ALWAYS, "#END SECOND INPUT BLOCK"); omniout_str(ALWAYS, "#BEGIN OVERRIDE BLOCK"); omniout_str(ALWAYS, "glob_h := 0.005 ;"); omniout_str(ALWAYS, "glob_display_interval := 0.1;"); omniout_str(ALWAYS, "glob_look_poles := true;"); omniout_str(ALWAYS, "glob_max_iter := 10000;"); omniout_str(ALWAYS, "glob_max_minutes := 10;"); omniout_str(ALWAYS, "#END OVERRIDE BLOCK"); omniout_str(ALWAYS, "!"); omniout_str(ALWAYS, "#BEGIN USER DEF BLOCK"); omniout_str(ALWAYS, "exact_soln_y := proc(x)"); omniout_str(ALWAYS, "return(1.0 / (x * x + 1.0));"); omniout_str(ALWAYS, "end;"); omniout_str(ALWAYS, ""); omniout_str(ALWAYS, ""); omniout_str(ALWAYS, "#END USER DEF BLOCK"); omniout_str(ALWAYS, "#######END OF ECHO OF PROBLEM#################"); glob_unchanged_h_cnt := 0; glob_warned := false; glob_warned2 := false; glob_small_float := 0.10*10^(-199); glob_smallish_float := 0.10*10^(-63); glob_large_float := 0.10*10^101; glob_almost_1 := 0.99; glob_log10_abserr := -8.0; glob_log10_relerr := -8.0; glob_hmax := 0.01; Digits := 32; max_terms := 30; glob_max_terms := max_terms; glob_html_log := true; array_y := Array(0 .. max_terms + 1, []); array_x := Array(0 .. max_terms + 1, []); array_norms := Array(0 .. max_terms + 1, []); array_y_init := Array(0 .. max_terms + 1, []); array_last_rel_error := Array(0 .. max_terms + 1, []); array_tmp0 := Array(0 .. max_terms + 1, []); array_tmp1 := Array(0 .. max_terms + 1, []); array_tmp2 := Array(0 .. max_terms + 1, []); array_tmp3 := Array(0 .. max_terms + 1, []); array_tmp4 := Array(0 .. max_terms + 1, []); array_tmp5 := Array(0 .. max_terms + 1, []); array_tmp6 := Array(0 .. max_terms + 1, []); array_tmp7 := Array(0 .. max_terms + 1, []); array_tmp8 := Array(0 .. max_terms + 1, []); array_tmp9 := Array(0 .. max_terms + 1, []); array_type_pole := Array(0 .. max_terms + 1, []); array_m1 := Array(0 .. max_terms + 1, []); array_fact_1 := Array(0 .. max_terms + 1, []); array_pole := Array(0 .. max_terms + 1, []); array_1st_rel_error := Array(0 .. max_terms + 1, []); array_y_higher := Array(0 .. 3, 0 .. max_terms + 1, []); array_real_pole := Array(0 .. 2, 0 .. 4, []); array_complex_pole := Array(0 .. 2, 0 .. 4, []); array_y_higher_work2 := Array(0 .. 3, 0 .. max_terms + 1, []); array_fact_2 := Array(0 .. max_terms + 1, 0 .. max_terms + 1, []); array_poles := Array(0 .. 2, 0 .. 4, []); array_y_set_initial := Array(0 .. 3, 0 .. max_terms + 1, []); array_y_higher_work := Array(0 .. 3, 0 .. max_terms + 1, []); term := 1; while term <= max_terms do array_y[term] := 0.; term := term + 1 end do ; term := 1; while term <= max_terms do array_x[term] := 0.; term := term + 1 end do ; term := 1; while term <= max_terms do array_norms[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_y_init[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_last_rel_error[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp0[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp1[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp2[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp3[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp4[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp5[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp6[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp7[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp8[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_tmp9[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_type_pole[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_m1[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_fact_1[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_pole[term] := 0.; term := term + 1 end do; term := 1; while term <= max_terms do array_1st_rel_error[term] := 0.; term := term + 1 end do; ord := 1; while ord <= 2 do term := 1; while term <= max_terms do array_y_higher[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 1 do term := 1; while term <= 3 do array_real_pole[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 1 do term := 1; while term <= 3 do array_complex_pole[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 2 do term := 1; while term <= max_terms do array_y_higher_work2[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= max_terms do term := 1; while term <= max_terms do array_fact_2[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 1 do term := 1; while term <= 3 do array_poles[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 2 do term := 1; while term <= max_terms do array_y_set_initial[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; ord := 1; while ord <= 2 do term := 1; while term <= max_terms do array_y_higher_work[ord, term] := 0.; term := term + 1 end do; ord := ord + 1 end do; array_x := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_x[term] := 0.; term := term + 1 end do; array_y := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_y[term] := 0.; term := term + 1 end do; array_tmp9 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp9[term] := 0.; term := term + 1 end do; array_tmp8 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp8[term] := 0.; term := term + 1 end do; array_tmp7 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp7[term] := 0.; term := term + 1 end do; array_tmp6 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp6[term] := 0.; term := term + 1 end do; array_tmp5 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp5[term] := 0.; term := term + 1 end do; array_tmp4 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp4[term] := 0.; term := term + 1 end do; array_tmp3 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp3[term] := 0.; term := term + 1 end do; array_tmp2 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp2[term] := 0.; term := term + 1 end do; array_tmp1 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp1[term] := 0.; term := term + 1 end do; array_tmp0 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_tmp0[term] := 0.; term := term + 1 end do; array_m1 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_m1[term] := 0.; term := term + 1 end do; array_const_1D0 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_const_1D0[term] := 0.; term := term + 1 end do; array_const_1D0[1] := 1.0; array_const_1 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_const_1[term] := 0.; term := term + 1 end do; array_const_1[1] := 1; array_const_2D0 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_const_2D0[term] := 0.; term := term + 1 end do; array_const_2D0[1] := 2.0; array_const_0D0 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms + 1 do array_const_0D0[term] := 0.; term := term + 1 end do; array_const_0D0[1] := 0.; array_m1 := Array(1 .. max_terms + 2, []); term := 1; while term <= max_terms do array_m1[term] := 0.; term := term + 1 end do; array_m1[1] := -1.0; iiif := 0; while iiif <= glob_max_terms do jjjf := 0; while jjjf <= glob_max_terms do array_fact_1[iiif] := 0; array_fact_2[iiif, jjjf] := 0; jjjf := jjjf + 1 end do; iiif := iiif + 1 end do; x_start := -2.0; x_end := 1.0; array_y_init[1] := exact_soln_y(x_start); glob_h := 0.1; glob_look_poles := true; glob_max_iter := 50; glob_h := 0.005; glob_display_interval := 0.1; glob_look_poles := true; glob_max_iter := 10000; glob_max_minutes := 10; 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 0. < glob_h then glob_neg_h := false; glob_display_interval := omniabs(glob_display_interval) else glob_neg_h := true; glob_display_interval := -omniabs(glob_display_interval) end if; chk_data(); 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; if glob_html_log then html_log_file := fopen("html/entry.html", WRITE, TEXT) end if; omniout_str(ALWAYS, "START of Soultion"); array_x[1] := x_start; array_x[2] := glob_h; glob_next_display := x_start; order_diff := 1; term_no := 1; while term_no <= order_diff do 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 do; rows := order_diff; r_order := 1; while r_order <= rows do term_no := 1; while term_no <= rows - r_order + 1 do 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 do; r_order := r_order + 1 end do; current_iter := 1; glob_clock_start_sec := elapsed_time_seconds(); if glob_small_float < omniabs(array_y_higher[1, 1]) then tmp := omniabs(array_y_higher[1, 1]); log10norm := log10(tmp); if log10norm < glob_log10normmin then glob_log10normmin := log10norm end if end if; 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 not_reached_end(array_x[1], x_end) and convfloat(glob_clock_sec) - convfloat(glob_orig_start_sec) < convfloat(glob_max_sec) do if reached_interval() then omniout_str(INFO, " "); omniout_str(INFO, "TOP MAIN SOLVE Loop") end if; glob_iter := glob_iter + 1; glob_clock_sec := elapsed_time_seconds(); glob_current_iter := glob_current_iter + 1; atomall(); if glob_look_poles then check_for_pole() end if; if reached_interval() then glob_next_display := glob_next_display + glob_display_interval end if; array_x[1] := array_x[1] + glob_h; array_x[2] := glob_h; order_diff := 1; ord := 2; calc_term := 1; iii := glob_max_terms; while calc_term <= iii do 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 do; temp_sum := 0.; ord := 2; calc_term := 1; iii := glob_max_terms; while calc_term <= iii do temp_sum := temp_sum + array_y_higher_work[ord, iii]; iii := iii - 1 end do; array_y_higher_work2[ord, calc_term] := temp_sum*expt(glob_h, calc_term - 1)/factorial_1(calc_term - 1) ; ord := 1; calc_term := 2; iii := glob_max_terms; while calc_term <= iii do 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 do; temp_sum := 0.; ord := 1; calc_term := 2; iii := glob_max_terms; while calc_term <= iii do temp_sum := temp_sum + array_y_higher_work[ord, iii]; iii := iii - 1 end do; array_y_higher_work2[ord, calc_term] := temp_sum*expt(glob_h, calc_term - 1)/factorial_1(calc_term - 1) ; ord := 1; calc_term := 1; iii := glob_max_terms; while calc_term <= iii do 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 do; temp_sum := 0.; ord := 1; calc_term := 1; iii := glob_max_terms; while calc_term <= iii do temp_sum := temp_sum + array_y_higher_work[ord, iii]; iii := iii - 1 end do; array_y_higher_work2[ord, calc_term] := temp_sum*expt(glob_h, calc_term - 1)/factorial_1(calc_term - 1) ; term_no := glob_max_terms; while 1 <= term_no do array_y[term_no] := array_y_higher_work2[1, term_no]; ord := 1; while ord <= order_diff do array_y_higher[ord, term_no] := array_y_higher_work2[ord, term_no]; ord := ord + 1 end do; term_no := term_no - 1 end do; display_alot(current_iter) end do; omniout_str(ALWAYS, "Finished!"); if glob_max_iter <= glob_iter then omniout_str(ALWAYS, "Maximum Iterations Reached before Solution Completed!") end if; if convfloat(glob_max_sec) <= elapsed_time_seconds() - convfloat(glob_orig_start_sec) then omniout_str(ALWAYS, "Maximum Time Reached before Solution Completed!") end if; glob_clock_sec := elapsed_time_seconds(); omniout_str(INFO, "diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0)\ /( x * x + 1.0);"); omniout_int(INFO, "Iterations ", 32, glob_iter, 4, " "); prog_report(x_start, x_end); if glob_html_log then logstart(html_log_file); logitem_str(html_log_file, "2012-09-21T03:44:08-05:00"); logitem_str(html_log_file, "Maple"); logitem_str(html_log_file, "sing4"); logitem_str(html_log_file, "diff ( y , x , 1 ) = m1 * 2.0 * x / (\ x * x + 1.0) /( x * x + 1.0);"); logitem_float(html_log_file, x_start); logitem_float(html_log_file, x_end); logitem_float(html_log_file, array_x[1]); logitem_float(html_log_file, glob_h); logitem_integer(html_log_file, Digits); 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 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 if; logitem_time(html_log_file, convfloat(glob_clock_sec)); if glob_percent_done < 100.0 then logitem_time(html_log_file, convfloat(glob_optimal_expect_sec)) ; 0 else logitem_str(html_log_file, "Done"); 0 end if; log_revs(html_log_file, " 130 "); logitem_str(html_log_file, "sing4 diffeq.mxt"); logitem_str(html_log_file, "sing4 maple results"); logitem_str(html_log_file, "c c++ Maple and Maxima"); logend(html_log_file) end if; if glob_html_log then fclose(html_log_file) end if end proc > main(); ##############ECHO OF PROBLEM################# ##############temp/sing4postode.ode################# diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0) /( x * x + 1.0); ! #BEGIN FIRST INPUT BLOCK Digits := 32; max_terms := 30; ! #END FIRST INPUT BLOCK #BEGIN SECOND INPUT BLOCK x_start := -2.0; x_end := 1.0; array_y_init[0 + 1] := exact_soln_y(x_start); glob_h := 0.1; glob_look_poles := true; glob_max_iter := 50; #END SECOND INPUT BLOCK #BEGIN OVERRIDE BLOCK glob_h := 0.005 ; glob_display_interval := 0.1; glob_look_poles := true; glob_max_iter := 10000; glob_max_minutes := 10; #END OVERRIDE BLOCK ! #BEGIN USER DEF BLOCK exact_soln_y := proc(x) return(1.0 / (x * x + 1.0)); end; #END USER DEF BLOCK #######END OF ECHO OF PROBLEM################# START of Soultion x[1] = -2 y[1] (analytic) = 0.2 y[1] (numeric) = 0.2 absolute error = 0 relative error = 0 % Correct digits = 32 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 2.259 Order of pole = 3.572 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=3.8MB, alloc=3.1MB, time=0.17 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.9 y[1] (analytic) = 0.21691973969631236442516268980477 y[1] (numeric) = 0.21691973969630443866369260757422 absolute error = 7.92576147008223055e-15 relative error = 3.6537760377079082835500000000000e-12 % Correct digits = 13 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 2.166 Order of pole = 3.531 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=7.6MB, alloc=4.3MB, time=0.38 Complex estimate of poles used x[1] = -1.8 y[1] (analytic) = 0.23584905660377358490566037735849 y[1] (numeric) = 0.23584905660374624729978184361677 absolute error = 2.733760587853374172e-14 relative error = 1.1591144892498306489280000000000e-11 % Correct digits = 12 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 2.076 Order of pole = 3.519 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.7 y[1] (analytic) = 0.25706940874035989717223650385604 y[1] (numeric) = 0.25706940874029419499109837166242 absolute error = 6.570218113813219362e-14 relative error = 2.5558148462733423318180000000000e-11 % Correct digits = 12 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.991 Order of pole = 3.559 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=11.4MB, alloc=4.4MB, time=0.59 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.6 y[1] (analytic) = 0.28089887640449438202247191011236 y[1] (numeric) = 0.28089887640436011716597361078652 absolute error = 1.3426485649829932584e-13 relative error = 4.7798288913394559999040000000000e-11 % Correct digits = 12 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.908 Order of pole = 3.609 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=15.2MB, alloc=4.5MB, time=0.80 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.5 y[1] (analytic) = 0.30769230769230769230769230769231 y[1] (numeric) = 0.30769230769205838556101444933018 absolute error = 2.4930674667785836213e-13 relative error = 8.1024692670303967692249999999999e-11 % Correct digits = 12 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.823 Order of pole = 3.599 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=19.0MB, alloc=4.5MB, time=1.01 Complex estimate of poles used Complex estimate of poles used x[1] = -1.4 y[1] (analytic) = 0.33783783783783783783783783783784 y[1] (numeric) = 0.33783783783740476463246066517914 absolute error = 4.3307320537717265870e-13 relative error = 1.2818966879164310697520000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.737 Order of pole = 3.555 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.3 y[1] (analytic) = 0.37174721189591078066914498141264 y[1] (numeric) = 0.37174721189519760188339366217278 absolute error = 7.1317878575131923986e-13 relative error = 1.9184509336710487552234000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.657 Order of pole = 3.577 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=22.8MB, alloc=4.6MB, time=1.22 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.2 y[1] (analytic) = 0.40983606557377049180327868852459 y[1] (numeric) = 0.40983606557265258953839574083057 absolute error = 1.11790226488294769402e-12 relative error = 2.7276815263143923734088000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.581 Order of pole = 3.622 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=26.7MB, alloc=4.6MB, time=1.44 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -1.1 y[1] (analytic) = 0.45248868778280542986425339366516 y[1] (numeric) = 0.45248868778114286762215356253655 absolute error = 1.66256224209983112861e-12 relative error = 3.6742625550406267942281000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.502 Order of pole = 3.581 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=30.5MB, alloc=4.6MB, time=1.65 Complex estimate of poles used Complex estimate of poles used x[1] = -1 y[1] (analytic) = 0.5 y[1] (numeric) = 0.49999999999768035292474613683826 absolute error = 2.31964707525386316174e-12 relative error = 4.6392941505077263234800000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop NO POLE Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.9 y[1] (analytic) = 0.55248618784530386740331491712707 y[1] (numeric) = 0.55248618784233847803421387618005 absolute error = 2.96538936910104094702e-12 relative error = 5.3673547580728841141062000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.361 Order of pole = 3.618 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=34.3MB, alloc=4.6MB, time=1.86 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.8 y[1] (analytic) = 0.60975609756097560975609756097561 y[1] (numeric) = 0.60975609755766977486349185281029 absolute error = 3.30583489260570816532e-12 relative error = 5.4215692238733613911248000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.294 Order of pole = 3.588 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=38.1MB, alloc=4.6MB, time=2.08 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.7 y[1] (analytic) = 0.67114093959731543624161073825503 y[1] (numeric) = 0.67114093959449529862674910383706 absolute error = 2.82013761486163441797e-12 relative error = 4.2020050461438352827753000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.235 Order of pole = 3.624 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=41.9MB, alloc=4.6MB, time=2.29 Complex estimate of poles used Complex estimate of poles used x[1] = -0.6 y[1] (analytic) = 0.73529411764705882352941176470588 y[1] (numeric) = 0.73529411764623195180800913640538 absolute error = 8.2687172140262830050e-13 relative error = 1.1245455411075744886800000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.179 Order of pole = 3.602 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.5 y[1] (analytic) = 0.8 y[1] (numeric) = 0.80000000000315826862117216729667 absolute error = 3.15826862117216729667e-12 relative error = 3.9478357764652091208375000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.131 Order of pole = 3.608 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=45.7MB, alloc=4.6MB, time=2.52 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.4 y[1] (analytic) = 0.86206896551724137931034482758621 y[1] (numeric) = 0.86206896552594169397457884789588 absolute error = 8.70031466423402030967e-12 relative error = 1.0092365010511463559217200000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.09 Order of pole = 3.631 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=49.5MB, alloc=4.6MB, time=2.73 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.3 y[1] (analytic) = 0.91743119266055045871559633027523 y[1] (numeric) = 0.91743119267429724476347706972234 absolute error = 1.374678604788073944711e-11 relative error = 1.4983996792190005997349900000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.056 Order of pole = 3.607 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=53.4MB, alloc=4.6MB, time=2.95 Complex estimate of poles used Complex estimate of poles used x[1] = -0.2 y[1] (analytic) = 0.96153846153846153846153846153846 y[1] (numeric) = 0.96153846155344780304495188290868 absolute error = 1.498626458341342137022e-11 relative error = 1.5585715166749958225028800000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.031 Order of pole = 3.602 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = -0.1 y[1] (analytic) = 0.99009900990099009900990099009901 y[1] (numeric) = 0.99009900991102274754480526476584 absolute error = 1.003264853490427466683e-11 relative error = 1.0132975020253317413498300000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.017 Order of pole = 3.622 x[1] = -0.095 y[1] (analytic) = 0.99105572210797552092366393300463 y[1] (numeric) = 0.99105572211760086075355954947605 absolute error = 9.62533982989561647142e-12 relative error = 9.7122085218604244100745655000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.016 Order of pole = 3.618 Complex estimate of poles used Complex estimate of poles used memory used=57.2MB, alloc=4.6MB, time=3.16 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=61.0MB, alloc=4.7MB, time=3.38 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.1 y[1] (analytic) = 0.99009900990099009900990099009901 y[1] (numeric) = 0.99009900989110081638448970357177 absolute error = 9.88928262541128652724e-12 relative error = 9.9881754516653993925124000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.017 Order of pole = 3.622 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=64.8MB, alloc=4.7MB, time=3.59 Complex estimate of poles used Complex estimate of poles used x[1] = 0.2 y[1] (analytic) = 0.96153846153846153846153846153846 y[1] (numeric) = 0.9615384615234269618043934255169 absolute error = 1.503457665714503602156e-11 relative error = 1.5635959723430837462422400000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.031 Order of pole = 3.602 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.3 y[1] (analytic) = 0.91743119266055045871559633027523 y[1] (numeric) = 0.91743119264658960514118009493665 absolute error = 1.396085357441623533858e-11 relative error = 1.5217330396113696519052200000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.056 Order of pole = 3.607 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=68.6MB, alloc=4.7MB, time=3.82 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.4 y[1] (analytic) = 0.86206896551724137931034482758621 y[1] (numeric) = 0.86206896550826261398421514428247 absolute error = 8.97876532612968330374e-12 relative error = 1.0415367778310432632338400000000e-09 % Correct digits = 10 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.09 Order of pole = 3.631 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=72.4MB, alloc=4.7MB, time=4.04 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.5 y[1] (analytic) = 0.8 y[1] (numeric) = 0.79999999999658745428432425478471 absolute error = 3.41254571567574521529e-12 relative error = 4.2656821445946815191125000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.131 Order of pole = 3.608 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=76.2MB, alloc=4.7MB, time=4.25 Complex estimate of poles used Complex estimate of poles used x[1] = 0.6 y[1] (analytic) = 0.73529411764705882352941176470588 y[1] (numeric) = 0.73529411764769205845029645783177 absolute error = 6.3323492088469312589e-13 relative error = 8.6119949240318265121040000000000e-11 % Correct digits = 12 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.179 Order of pole = 3.602 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.7 y[1] (analytic) = 0.67114093959731543624161073825503 y[1] (numeric) = 0.67114093959999790236948443820764 absolute error = 2.68246612787369995261e-12 relative error = 3.9968745305318129293889000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.235 Order of pole = 3.624 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=80.1MB, alloc=4.7MB, time=4.47 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.8 y[1] (analytic) = 0.60975609756097560975609756097561 y[1] (numeric) = 0.60975609756417929161267830169753 absolute error = 3.20368185658074072192e-12 relative error = 5.2540382447924147839488000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.294 Order of pole = 3.588 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=83.9MB, alloc=4.7MB, time=4.69 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used x[1] = 0.9 y[1] (analytic) = 0.55248618784530386740331491712707 y[1] (numeric) = 0.55248618784818336509487050630935 absolute error = 2.87949769155558918228e-12 relative error = 5.2118908217156164199268000000000e-10 % Correct digits = 11 h = 0.005 TOP MAIN SOLVE Loop Complex estimate of poles used Radius of convergence = 1.361 Order of pole = 3.618 Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used Complex estimate of poles used memory used=87.7MB, alloc=4.7MB, time=4.90 Complex estimate of poles used Complex estimate of poles used x[1] = 1 y[1] (analytic) = 0.5 y[1] (numeric) = 0.50000000000223764459231002030033 absolute error = 2.23764459231002030033e-12 relative error = 4.4752891846200406006600000000000e-10 % Correct digits = 11 h = 0.005 Finished! diff ( y , x , 1 ) = m1 * 2.0 * x / (x * x + 1.0) /( x * x + 1.0); Iterations = 600 Total Elapsed Time = 4 Seconds Elapsed Time(since restart) = 4 Seconds Time to Timeout = 9 Minutes 55 Seconds Percent Done = 100.2 % > quit memory used=88.0MB, alloc=4.7MB, time=4.92