|\^/| 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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> 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 CONST - LINEAR $eq_no = 1 i = 1
> array_tmp1[1] := array_const_0D1[1] * array_x[1];
> #emit pre add LINEAR - CONST $eq_no = 1 i = 1
> array_tmp2[1] := array_tmp1[1] + array_const_0D2[1];
> #emit pre mult CONST - LINEAR $eq_no = 1 i = 1
> array_tmp3[1] := array_const_0D2[1] * array_x[1];
> #emit pre add LINEAR - CONST $eq_no = 1 i = 1
> array_tmp4[1] := array_tmp3[1] + array_const_0D3[1];
> #emit pre div LINEAR - LINEAR $eq_no = 1 i = 1
> array_tmp5[1] := array_tmp2[1] / array_tmp4[1];
> #emit pre add CONST FULL $eq_no = 1 i = 1
> array_tmp6[1] := array_const_0D0[1] + array_tmp5[1];
> #emit pre assign xxx $eq_no = 1 i = 1 $min_hdrs = 5
> if ( not array_y_set_initial[1,2]) then # if number 1
> if (1 <= glob_max_terms) then # if number 2
> temporary := array_tmp6[1] * expt(glob_h , (1)) * factorial_3(0,1);
> array_y[2] := temporary;
> array_y_higher[1,2] := temporary;
> temporary := temporary / glob_h * (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 CONST - LINEAR $eq_no = 1 i = 2
> array_tmp1[2] := array_const_0D1[1] * array_x[2];
> #emit pre add LINEAR - CONST $eq_no = 1 i = 2
> array_tmp2[2] := array_tmp1[2];
> #emit pre mult CONST - LINEAR $eq_no = 1 i = 2
> array_tmp3[2] := array_const_0D2[1] * array_x[2];
> #emit pre add LINEAR - CONST $eq_no = 1 i = 2
> array_tmp4[2] := array_tmp3[2];
> #emit pre div LINEAR - LINEAR $eq_no = 1 i = 2
> array_tmp5[2] := (array_tmp2[2] - array_tmp5[1] * array_tmp4[2]) / array_tmp4[1];
> #emit pre add CONST FULL $eq_no = 1 i = 2
> array_tmp6[2] := array_tmp5[2];
> #emit pre assign xxx $eq_no = 1 i = 2 $min_hdrs = 5
> if ( not array_y_set_initial[1,3]) then # if number 1
> if (2 <= glob_max_terms) then # if number 2
> temporary := array_tmp6[2] * expt(glob_h , (1)) * factorial_3(1,2);
> array_y[3] := temporary;
> array_y_higher[1,3] := temporary;
> temporary := temporary / glob_h * (2.0);
> array_y_higher[2,2] := temporary
> ;
> fi;# end if 2
> ;
> fi;# end if 1
> ;
> kkk := 3;
> #END ATOMHDR2
> #BEGIN ATOMHDR3
> #emit pre div LINEAR - LINEAR $eq_no = 1 i = 3
> array_tmp5[3] := - array_tmp5[2] * array_tmp4[2] / array_tmp4[1];
> #emit pre add CONST FULL $eq_no = 1 i = 3
> array_tmp6[3] := array_tmp5[3];
> #emit pre assign xxx $eq_no = 1 i = 3 $min_hdrs = 5
> if ( not array_y_set_initial[1,4]) then # if number 1
> if (3 <= glob_max_terms) then # if number 2
> temporary := array_tmp6[3] * expt(glob_h , (1)) * factorial_3(2,3);
> array_y[4] := temporary;
> array_y_higher[1,4] := temporary;
> temporary := temporary / glob_h * (2.0);
> array_y_higher[2,3] := temporary
> ;
> fi;# end if 2
> ;
> fi;# end if 1
> ;
> kkk := 4;
> #END ATOMHDR3
> #BEGIN ATOMHDR4
> #emit pre div LINEAR - LINEAR $eq_no = 1 i = 4
> array_tmp5[4] := - array_tmp5[3] * array_tmp4[2] / array_tmp4[1];
> #emit pre add CONST FULL $eq_no = 1 i = 4
> array_tmp6[4] := array_tmp5[4];
> #emit pre assign xxx $eq_no = 1 i = 4 $min_hdrs = 5
> if ( not array_y_set_initial[1,5]) then # if number 1
> if (4 <= glob_max_terms) then # if number 2
> temporary := array_tmp6[4] * expt(glob_h , (1)) * factorial_3(3,4);
> array_y[5] := temporary;
> array_y_higher[1,5] := temporary;
> temporary := temporary / glob_h * (2.0);
> array_y_higher[2,4] := temporary
> ;
> fi;# end if 2
> ;
> fi;# end if 1
> ;
> kkk := 5;
> #END ATOMHDR4
> #BEGIN ATOMHDR5
> #emit pre div LINEAR - LINEAR $eq_no = 1 i = 5
> array_tmp5[5] := - array_tmp5[4] * array_tmp4[2] / array_tmp4[1];
> #emit pre add CONST FULL $eq_no = 1 i = 5
> array_tmp6[5] := array_tmp5[5];
> #emit pre assign xxx $eq_no = 1 i = 5 $min_hdrs = 5
> if ( not array_y_set_initial[1,6]) then # if number 1
> if (5 <= glob_max_terms) then # if number 2
> temporary := array_tmp6[5] * expt(glob_h , (1)) * factorial_3(4,5);
> array_y[6] := temporary;
> array_y_higher[1,6] := temporary;
> temporary := temporary / glob_h * (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 div LINEAR - LINEAR (NOP) $eq_no = 1 i = 1
> array_tmp5[kkk] := - array_tmp5[kkk-1] * array_tmp4[2] / array_tmp4[1];
> #emit NOT FULL - FULL add $eq_no = 1
> array_tmp6[kkk] := array_tmp5[kkk];
> #emit assign $eq_no = 1
> order_d := 1;
> if (kkk + order_d + 1 <= glob_max_terms) then # if number 1
> if ( not array_y_set_initial[1,kkk + order_d]) then # if number 2
> temporary := array_tmp6[kkk] * expt(glob_h , (order_d)) / factorial_3((kkk - 1),(kkk + order_d - 1));
> array_y[kkk + order_d] := temporary;
> array_y_higher[1,kkk + order_d] := temporary;
> term := kkk + order_d - 1;
> adj2 := 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
array_y_set_initial, array_y_higher_work, glob_last;
array_tmp1[1] := array_const_0D1[1]*array_x[1];
array_tmp2[1] := array_tmp1[1] + array_const_0D2[1];
array_tmp3[1] := array_const_0D2[1]*array_x[1];
array_tmp4[1] := array_tmp3[1] + array_const_0D3[1];
array_tmp5[1] := array_tmp2[1]/array_tmp4[1];
array_tmp6[1] := array_const_0D0[1] + array_tmp5[1];
if not array_y_set_initial[1, 2] then
if 1 <= glob_max_terms then
temporary := array_tmp6[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_const_0D1[1]*array_x[2];
array_tmp2[2] := array_tmp1[2];
array_tmp3[2] := array_const_0D2[1]*array_x[2];
array_tmp4[2] := array_tmp3[2];
array_tmp5[2] :=
(array_tmp2[2] - array_tmp5[1]*array_tmp4[2])/array_tmp4[1];
array_tmp6[2] := array_tmp5[2];
if not array_y_set_initial[1, 3] then
if 2 <= glob_max_terms then
temporary := array_tmp6[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_tmp5[3] := -array_tmp5[2]*array_tmp4[2]/array_tmp4[1];
array_tmp6[3] := array_tmp5[3];
if not array_y_set_initial[1, 4] then
if 3 <= glob_max_terms then
temporary := array_tmp6[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_tmp5[4] := -array_tmp5[3]*array_tmp4[2]/array_tmp4[1];
array_tmp6[4] := array_tmp5[4];
if not array_y_set_initial[1, 5] then
if 4 <= glob_max_terms then
temporary := array_tmp6[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_tmp5[5] := -array_tmp5[4]*array_tmp4[2]/array_tmp4[1];
array_tmp6[5] := array_tmp5[5];
if not array_y_set_initial[1, 6] then
if 5 <= glob_max_terms then
temporary := array_tmp6[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_tmp5[kkk] := -array_tmp5[kkk - 1]*array_tmp4[2]/array_tmp4[1]
;
array_tmp6[kkk] := array_tmp5[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_tmp6[kkk]*expt(glob_h, order_d)/
factorial_3(kkk - 1, kkk + order_d - 1);
array_y[kkk + order_d] := temporary;
array_y_higher[1, kkk + order_d] := temporary;
term := kkk + order_d - 1;
adj2 := 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(0.5 * x + 0.25 * ln(2.0 * x + 3.0));
> end;
exact_soln_y := proc(x) return 0.5*x + 0.25*ln(2.0*x + 3.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
> DEBUGMASSIVE,
> INFO,
> glob_max_terms,
> DEBUGL,
> ALWAYS,
> glob_iolevel,
> #Top Generate Globals Decl
> MAX_UNCHANGED,
> glob_warned,
> glob_display_interval,
> glob_disp_incr,
> glob_reached_optimal_h,
> sec_in_minute,
> glob_display_flag,
> glob_smallish_float,
> glob_max_trunc_err,
> glob_log10_abserr,
> hours_in_day,
> glob_normmax,
> glob_iter,
> glob_unchanged_h_cnt,
> glob_max_rel_trunc_err,
> glob_neg_h,
> glob_hmin,
> glob_clock_sec,
> centuries_in_millinium,
> days_in_year,
> djd_debug,
> glob_max_iter,
> glob_last_good_h,
> glob_dump,
> glob_subiter_method,
> glob_abserr,
> glob_max_minutes,
> glob_max_sec,
> glob_max_hours,
> glob_dump_analytic,
> glob_optimal_done,
> glob_not_yet_start_msg,
> glob_not_yet_finished,
> years_in_century,
> glob_good_digits,
> glob_optimal_expect_sec,
> glob_log10abserr,
> glob_start,
> glob_optimal_clock_start_sec,
> glob_relerr,
> glob_next_display,
> glob_look_poles,
> glob_large_float,
> min_in_hour,
> glob_html_log,
> glob_log10relerr,
> glob_current_iter,
> glob_curr_iter_when_opt,
> glob_log10_relerr,
> glob_hmin_init,
> glob_hmax,
> djd_debug2,
> glob_log10normmin,
> glob_small_float,
> glob_optimal_start,
> glob_h,
> glob_percent_done,
> glob_no_eqs,
> glob_orig_start_sec,
> glob_warned2,
> glob_initial_pass,
> glob_clock_start_sec,
> glob_almost_1,
> glob_max_opt_iter,
> #Bottom Generate Globals Decl
> #BEGIN CONST
> array_const_1,
> array_const_0D3,
> array_const_0D2,
> array_const_0D1,
> array_const_0D0,
> #END CONST
> array_type_pole,
> array_m1,
> array_1st_rel_error,
> array_pole,
> array_last_rel_error,
> array_norms,
> array_y_init,
> array_y,
> array_x,
> array_tmp0,
> array_tmp1,
> array_tmp2,
> array_tmp3,
> array_tmp4,
> array_tmp5,
> array_tmp6,
> array_fact_1,
> array_y_higher_work2,
> array_poles,
> array_real_pole,
> array_complex_pole,
> array_y_higher,
> array_fact_2,
> array_y_set_initial,
> array_y_higher_work,
> glob_last;
> glob_last;
> ALWAYS := 1;
> INFO := 2;
> DEBUGL := 3;
> DEBUGMASSIVE := 4;
> glob_iolevel := INFO;
> DEBUGMASSIVE := 4;
> INFO := 2;
> glob_max_terms := 30;
> DEBUGL := 3;
> ALWAYS := 1;
> glob_iolevel := 5;
> MAX_UNCHANGED := 10;
> glob_warned := false;
> glob_display_interval := 0.0;
> glob_disp_incr := 0.1;
> glob_reached_optimal_h := false;
> sec_in_minute := 60;
> glob_display_flag := true;
> glob_smallish_float := 0.1e-100;
> glob_max_trunc_err := 0.1e-10;
> glob_log10_abserr := 0.1e-10;
> hours_in_day := 24;
> glob_normmax := 0.0;
> glob_iter := 0;
> glob_unchanged_h_cnt := 0;
> glob_max_rel_trunc_err := 0.1e-10;
> glob_neg_h := false;
> glob_hmin := 0.00000000001;
> glob_clock_sec := 0.0;
> centuries_in_millinium := 10;
> days_in_year := 365;
> djd_debug := true;
> glob_max_iter := 1000;
> glob_last_good_h := 0.1;
> glob_dump := false;
> glob_subiter_method := 3;
> glob_abserr := 0.1e-10;
> glob_max_minutes := 0.0;
> glob_max_sec := 10000.0;
> glob_max_hours := 0.0;
> glob_dump_analytic := false;
> glob_optimal_done := false;
> glob_not_yet_start_msg := true;
> glob_not_yet_finished := true;
> years_in_century := 100;
> glob_good_digits := 0;
> glob_optimal_expect_sec := 0.1;
> glob_log10abserr := 0.0;
> glob_start := 0;
> glob_optimal_clock_start_sec := 0.0;
> glob_relerr := 0.1e-10;
> glob_next_display := 0.0;
> glob_look_poles := false;
> glob_large_float := 9.0e100;
> min_in_hour := 60;
> glob_html_log := true;
> glob_log10relerr := 0.0;
> glob_current_iter := 0;
> glob_curr_iter_when_opt := 0;
> glob_log10_relerr := 0.1e-10;
> glob_hmin_init := 0.001;
> glob_hmax := 1.0;
> djd_debug2 := true;
> glob_log10normmin := 0.1;
> glob_small_float := 0.1e-50;
> glob_optimal_start := 0.0;
> glob_h := 0.1;
> glob_percent_done := 0.0;
> glob_no_eqs := 0;
> glob_orig_start_sec := 0.0;
> glob_warned2 := false;
> glob_initial_pass := true;
> glob_clock_start_sec := 0.0;
> glob_almost_1 := 0.9990;
> glob_max_opt_iter := 10;
> #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/div_lin_linpostode.ode#################");
> omniout_str(ALWAYS,"diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);");
> 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 := 0.1;");
> omniout_str(ALWAYS,"x_end := 5.0 ;");
> omniout_str(ALWAYS,"array_y_init[0 + 1] := exact_soln_y(x_start);");
> omniout_str(ALWAYS,"glob_h := 0.05;");
> omniout_str(ALWAYS,"glob_look_poles := true;");
> omniout_str(ALWAYS,"glob_max_iter := 1000000;");
> omniout_str(ALWAYS,"#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(0.5 * x + 0.25 * ln(2.0 * x + 3.0));");
> omniout_str(ALWAYS,"end;");
> omniout_str(ALWAYS,"#END USER DEF BLOCK");
> omniout_str(ALWAYS,"#######END OF ECHO OF PROBLEM#################");
> glob_unchanged_h_cnt := 0;
> glob_warned := false;
> glob_warned2 := false;
> glob_small_float := 1.0e-200;
> glob_smallish_float := 1.0e-64;
> glob_large_float := 1.0e100;
> glob_almost_1 := 0.99;
> glob_log10_abserr := -8.0;
> glob_log10_relerr := -8.0;
> glob_hmax := 0.01;
> #BEGIN FIRST INPUT BLOCK
> #BEGIN FIRST INPUT BLOCK
> Digits := 32;
> 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_type_pole:= Array(0..(max_terms + 1),[]);
> array_m1:= Array(0..(max_terms + 1),[]);
> array_1st_rel_error:= Array(0..(max_terms + 1),[]);
> array_pole:= Array(0..(max_terms + 1),[]);
> array_last_rel_error:= Array(0..(max_terms + 1),[]);
> array_norms:= Array(0..(max_terms + 1),[]);
> array_y_init:= Array(0..(max_terms + 1),[]);
> array_y:= Array(0..(max_terms + 1),[]);
> array_x:= 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_fact_1:= Array(0..(max_terms + 1),[]);
> array_y_higher_work2 := Array(0..(2+ 1) ,(0..max_terms+ 1),[]);
> array_poles := Array(0..(1+ 1) ,(0..3+ 1),[]);
> array_real_pole := Array(0..(1+ 1) ,(0..3+ 1),[]);
> array_complex_pole := Array(0..(1+ 1) ,(0..3+ 1),[]);
> array_y_higher := Array(0..(2+ 1) ,(0..max_terms+ 1),[]);
> array_fact_2 := Array(0..(max_terms+ 1) ,(0..max_terms+ 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_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_1st_rel_error[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_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_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_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_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_fact_1[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_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 <=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 <=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[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 <=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_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_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_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_0D3 := Array(1..(max_terms+1 + 1),[]);
> term := 1;
> while (term <= max_terms + 1) do # do number 2
> array_const_0D3[term] := 0.0;
> term := term + 1;
> od;# end do number 2
> ;
> array_const_0D3[1] := 0.3;
> array_const_0D2 := Array(1..(max_terms+1 + 1),[]);
> term := 1;
> while (term <= max_terms + 1) do # do number 2
> array_const_0D2[term] := 0.0;
> term := term + 1;
> od;# end do number 2
> ;
> array_const_0D2[1] := 0.2;
> array_const_0D1 := Array(1..(max_terms+1 + 1),[]);
> term := 1;
> while (term <= max_terms + 1) do # do number 2
> array_const_0D1[term] := 0.0;
> term := term + 1;
> od;# end do number 2
> ;
> array_const_0D1[1] := 0.1;
> 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 := 0.1;
> x_end := 5.0 ;
> array_y_init[0 + 1] := exact_soln_y(x_start);
> glob_h := 0.05;
> glob_look_poles := true;
> glob_max_iter := 1000000;
> #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 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);");
> 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-20T22:32:50-05:00")
> ;
> logitem_str(html_log_file,"Maple")
> ;
> logitem_str(html_log_file,"div_lin_lin")
> ;
> logitem_str(html_log_file,"diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);")
> ;
> 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,"div_lin_lin diffeq.mxt")
> ;
> logitem_str(html_log_file,"div_lin_lin 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 DEBUGMASSIVE, INFO, glob_max_terms, DEBUGL, ALWAYS, glob_iolevel,
MAX_UNCHANGED, glob_warned, glob_display_interval, glob_disp_incr,
glob_reached_optimal_h, sec_in_minute, glob_display_flag,
glob_smallish_float, glob_max_trunc_err, glob_log10_abserr, hours_in_day,
glob_normmax, glob_iter, glob_unchanged_h_cnt, glob_max_rel_trunc_err,
glob_neg_h, glob_hmin, glob_clock_sec, centuries_in_millinium, days_in_year,
djd_debug, glob_max_iter, glob_last_good_h, glob_dump, glob_subiter_method,
glob_abserr, glob_max_minutes, glob_max_sec, glob_max_hours,
glob_dump_analytic, glob_optimal_done, glob_not_yet_start_msg,
glob_not_yet_finished, years_in_century, glob_good_digits,
glob_optimal_expect_sec, glob_log10abserr, glob_start,
glob_optimal_clock_start_sec, glob_relerr, glob_next_display,
glob_look_poles, glob_large_float, min_in_hour, glob_html_log,
glob_log10relerr, glob_current_iter, glob_curr_iter_when_opt,
glob_log10_relerr, glob_hmin_init, glob_hmax, djd_debug2, glob_log10normmin,
glob_small_float, glob_optimal_start, glob_h, glob_percent_done,
glob_no_eqs, glob_orig_start_sec, glob_warned2, glob_initial_pass,
glob_clock_start_sec, glob_almost_1, glob_max_opt_iter, array_const_1,
array_const_0D3, array_const_0D2, array_const_0D1, array_const_0D0,
array_type_pole, array_m1, array_1st_rel_error, array_pole,
array_last_rel_error, array_norms, array_y_init, array_y, array_x,
array_tmp0, array_tmp1, array_tmp2, array_tmp3, array_tmp4, array_tmp5,
array_tmp6, array_fact_1, array_y_higher_work2, array_poles,
array_real_pole, array_complex_pole, array_y_higher, array_fact_2,
array_y_set_initial, array_y_higher_work, glob_last;
glob_last;
ALWAYS := 1;
INFO := 2;
DEBUGL := 3;
DEBUGMASSIVE := 4;
glob_iolevel := INFO;
DEBUGMASSIVE := 4;
INFO := 2;
glob_max_terms := 30;
DEBUGL := 3;
ALWAYS := 1;
glob_iolevel := 5;
MAX_UNCHANGED := 10;
glob_warned := false;
glob_display_interval := 0.;
glob_disp_incr := 0.1;
glob_reached_optimal_h := false;
sec_in_minute := 60;
glob_display_flag := true;
glob_smallish_float := 0.1*10^(-100);
glob_max_trunc_err := 0.1*10^(-10);
glob_log10_abserr := 0.1*10^(-10);
hours_in_day := 24;
glob_normmax := 0.;
glob_iter := 0;
glob_unchanged_h_cnt := 0;
glob_max_rel_trunc_err := 0.1*10^(-10);
glob_neg_h := false;
glob_hmin := 0.1*10^(-10);
glob_clock_sec := 0.;
centuries_in_millinium := 10;
days_in_year := 365;
djd_debug := true;
glob_max_iter := 1000;
glob_last_good_h := 0.1;
glob_dump := false;
glob_subiter_method := 3;
glob_abserr := 0.1*10^(-10);
glob_max_minutes := 0.;
glob_max_sec := 10000.0;
glob_max_hours := 0.;
glob_dump_analytic := false;
glob_optimal_done := false;
glob_not_yet_start_msg := true;
glob_not_yet_finished := true;
years_in_century := 100;
glob_good_digits := 0;
glob_optimal_expect_sec := 0.1;
glob_log10abserr := 0.;
glob_start := 0;
glob_optimal_clock_start_sec := 0.;
glob_relerr := 0.1*10^(-10);
glob_next_display := 0.;
glob_look_poles := false;
glob_large_float := 0.90*10^101;
min_in_hour := 60;
glob_html_log := true;
glob_log10relerr := 0.;
glob_current_iter := 0;
glob_curr_iter_when_opt := 0;
glob_log10_relerr := 0.1*10^(-10);
glob_hmin_init := 0.001;
glob_hmax := 1.0;
djd_debug2 := true;
glob_log10normmin := 0.1;
glob_small_float := 0.1*10^(-50);
glob_optimal_start := 0.;
glob_h := 0.1;
glob_percent_done := 0.;
glob_no_eqs := 0;
glob_orig_start_sec := 0.;
glob_warned2 := false;
glob_initial_pass := true;
glob_clock_start_sec := 0.;
glob_almost_1 := 0.9990;
glob_max_opt_iter := 10;
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/div_lin_linpostode.ode#################");
omniout_str(ALWAYS,
"diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);");
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 := 0.1;");
omniout_str(ALWAYS, "x_end := 5.0 ;");
omniout_str(ALWAYS, "array_y_init[0 + 1] := exact_soln_y(x_start);");
omniout_str(ALWAYS, "glob_h := 0.05;");
omniout_str(ALWAYS, "glob_look_poles := true;");
omniout_str(ALWAYS, "glob_max_iter := 1000000;");
omniout_str(ALWAYS, "#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(0.5 * x + 0.25 * ln(2.0 * x + 3.0));");
omniout_str(ALWAYS, "end;");
omniout_str(ALWAYS, "#END USER DEF BLOCK");
omniout_str(ALWAYS, "#######END OF ECHO OF PROBLEM#################");
glob_unchanged_h_cnt := 0;
glob_warned := false;
glob_warned2 := false;
glob_small_float := 0.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_type_pole := Array(0 .. max_terms + 1, []);
array_m1 := Array(0 .. max_terms + 1, []);
array_1st_rel_error := Array(0 .. max_terms + 1, []);
array_pole := Array(0 .. max_terms + 1, []);
array_last_rel_error := Array(0 .. max_terms + 1, []);
array_norms := Array(0 .. max_terms + 1, []);
array_y_init := Array(0 .. max_terms + 1, []);
array_y := Array(0 .. max_terms + 1, []);
array_x := 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_fact_1 := Array(0 .. max_terms + 1, []);
array_y_higher_work2 := Array(0 .. 3, 0 .. max_terms + 1, []);
array_poles := Array(0 .. 2, 0 .. 4, []);
array_real_pole := Array(0 .. 2, 0 .. 4, []);
array_complex_pole := Array(0 .. 2, 0 .. 4, []);
array_y_higher := Array(0 .. 3, 0 .. max_terms + 1, []);
array_fact_2 := Array(0 .. max_terms + 1, 0 .. max_terms + 1, []);
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_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_1st_rel_error[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_last_rel_error[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_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_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_fact_1[term] := 0.; term := term + 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 <= 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 <= 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[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 <= 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_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_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_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_0D3 := Array(1 .. max_terms + 2, []);
term := 1;
while term <= max_terms + 1 do
array_const_0D3[term] := 0.; term := term + 1
end do;
array_const_0D3[1] := 0.3;
array_const_0D2 := Array(1 .. max_terms + 2, []);
term := 1;
while term <= max_terms + 1 do
array_const_0D2[term] := 0.; term := term + 1
end do;
array_const_0D2[1] := 0.2;
array_const_0D1 := Array(1 .. max_terms + 2, []);
term := 1;
while term <= max_terms + 1 do
array_const_0D1[term] := 0.; term := term + 1
end do;
array_const_0D1[1] := 0.1;
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 := 0.1;
x_end := 5.0;
array_y_init[1] := exact_soln_y(x_start);
glob_h := 0.05;
glob_look_poles := true;
glob_max_iter := 1000000;
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 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);");
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-20T22:32:50-05:00");
logitem_str(html_log_file, "Maple");
logitem_str(html_log_file,
"div_lin_lin");
logitem_str(html_log_file,
"diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);");
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, "div_lin_lin diffeq.mxt");
logitem_str(html_log_file, "div_lin_lin 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/div_lin_linpostode.ode#################
diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);
!
#BEGIN FIRST INPUT BLOCK
Digits := 32;
max_terms := 30;
!
#END FIRST INPUT BLOCK
#BEGIN SECOND INPUT BLOCK
x_start := 0.1;
x_end := 5.0 ;
array_y_init[0 + 1] := exact_soln_y(x_start);
glob_h := 0.05;
glob_look_poles := true;
glob_max_iter := 1000000;
#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(0.5 * x + 0.25 * ln(2.0 * x + 3.0));
end;
#END USER DEF BLOCK
#######END OF ECHO OF PROBLEM#################
START of Soultion
x[1] = 0.1
y[1] (analytic) = 0.34078770245142021576704228815162
y[1] (numeric) = 0.34078770245142021576704228815162
absolute error = 0
relative error = 0 %
Correct digits = 32
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
x[1] = 0.2
y[1] (analytic) = 0.40594385790552892641219382116172
y[1] (numeric) = 0.40594385790550607055312690912536
absolute error = 2.285585906691203636e-14
relative error = 5.6303005013641668868058997285426e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 4.436
Order of pole = 3.905
Complex estimate of poles used
Complex estimate of poles used
memory used=3.8MB, alloc=3.1MB, time=0.33
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 0.3
y[1] (analytic) = 0.47023346136551607940174081551925
y[1] (numeric) = 0.47023346136547718276866985676055
absolute error = 3.889663307095875870e-14
relative error = 8.2717705707302073964404557138119e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
x[1] = 0.4
y[1] (analytic) = 0.53375026668308502135206702466542
y[1] (numeric) = 0.53375026668303464213907272455692
absolute error = 5.037921299430010850e-14
relative error = 9.4387237138773811128551414277264e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 6.186
Order of pole = 5.8
Complex estimate of poles used
memory used=7.6MB, alloc=4.3MB, time=0.71
Complex 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.5965735902799726547086160607291
y[1] (numeric) = 0.59657359027991390925845476073546
absolute error = 5.874545016129999364e-14
relative error = 9.8471422668460211291916226002373e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 1.825
Order of pole = 0.9163
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 0.6
y[1] (analytic) = 0.65877113132233065547495966178488
y[1] (numeric) = 0.65877113132226571668196779552373
absolute error = 6.493879299186626115e-14
relative error = 9.8575650790156296355582399805880e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=11.4MB, alloc=4.4MB, time=1.09
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 0.7
y[1] (analytic) = 0.72040113523105386971960409154928
y[1] (numeric) = 0.72040113523098427974642376540153
absolute error = 6.958997318032614775e-14
relative error = 9.6598922151901652648090758252908e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 1.31
Order of pole = 0.6772
Complex 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.781514075873762329051498374646
y[1] (numeric) = 0.78151407587368920013263856695076
absolute error = 7.312891885980769524e-14
relative error = 9.3573386734009614491018832139411e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
memory used=15.2MB, alloc=4.4MB, time=1.48
Complex estimate of poles used
Complex 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.84215397947846131126154556701772
y[1] (numeric) = 0.84215397947838545737245796361446
absolute error = 7.585388908760340326e-14
relative error = 9.0071282611024484959323545736935e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 2.085
Order of pole = 0.8773
Complex estimate of poles used
Complex 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
y[1] (analytic) = 0.90235947810852509365018983330655
y[1] (numeric) = 0.90235947810844711850309753201449
absolute error = 7.797514709230129206e-14
relative error = 8.6412509630583573656210864105732e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=19.0MB, alloc=4.4MB, time=1.87
Complex estimate of poles used
Complex estimate of poles used
x[1] = 1.1
y[1] (analytic) = 0.96216465639684541771749005744932
y[1] (numeric) = 0.96216465639676577460096412510352
absolute error = 7.964311652593234580e-14
relative error = 8.2774934618969720599045817084524e-12 %
Correct digits = 13
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
x[1] = 1.2
y[1] (analytic) = 1.0215997383925571748962440943854
y[1] (numeric) = 1.0215997383924762080213822941078
absolute error = 8.09668748618002776e-14
relative error = 7.9254987857767210459209000985628e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=22.8MB, alloc=4.5MB, time=2.27
Complex estimate of poles used
x[1] = 1.3
y[1] (analytic) = 1.0806916494352758873347644132833
y[1] (numeric) = 1.0806916494351938608253472841366
absolute error = 8.20265094171291467e-14
relative error = 7.5901862904180631800201073713229e-12 %
Correct digits = 13
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
x[1] = 1.4
y[1] (analytic) = 1.1394644793880934131456281747839
y[1] (numeric) = 1.139464479388010531641844471012
absolute error = 8.28815037837037719e-14
relative error = 7.2737242172052760692198742134341e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=26.7MB, alloc=4.5MB, time=2.68
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 1.5
y[1] (analytic) = 1.1979398673070137502031193395952
y[1] (numeric) = 1.197939867306930173689717667983
absolute error = 8.35765134016716122e-14
relative error = 6.9766868673928375244972727828884e-12 %
Correct digits = 13
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
x[1] = 1.6
y[1] (analytic) = 1.256137323012761467832101247829
y[1] (numeric) = 1.2561373230126773224349875592397
absolute error = 8.41453971136885893e-14
relative error = 6.6987418948647648381299903023852e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
memory used=30.5MB, alloc=4.5MB, time=3.11
Complex estimate of poles used
x[1] = 1.7
y[1] (analytic) = 1.3140744975914065431213503185162
y[1] (numeric) = 1.314074497591321929054493828926
absolute error = 8.46140668564895902e-14
relative error = 6.4390616370365917209447877635919e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 8.863
Order of pole = 4.335
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 1.8
y[1] (analytic) = 1.3717674122580949652141073704154
y[1] (numeric) = 1.3717674122580099626884571382857
absolute error = 8.50025256502321297e-14
relative error = 6.1965698332421892337843130367294e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
memory used=34.3MB, alloc=4.5MB, time=3.54
Complex estimate of poles used
Complex estimate of poles used
x[1] = 1.9
y[1] (analytic) = 1.4292306530455152537665018515263
y[1] (numeric) = 1.4292306530454299274253976070039
absolute error = 8.53263411042445224e-14
relative error = 5.9700889371792126872440531151096e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 2.491
Order of pole = 0.7683
Complex 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] = 2
y[1] (analytic) = 1.4864775372638283262763381858608
y[1] (numeric) = 1.4864775372637427285545882035145
absolute error = 8.55977217499823463e-14
relative error = 5.7584268584066726532456543238030e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 2.786
Order of pole = 0.8169
Complex estimate of poles used
memory used=38.1MB, alloc=4.5MB, time=3.95
Complex estimate of poles used
x[1] = 2.1
y[1] (analytic) = 1.5435202565055024067560488458838
y[1] (numeric) = 1.5435202565054165804452039754252
absolute error = 8.58263108448704586e-14
relative error = 5.5604265951895851462938770610419e-12 %
Correct digits = 13
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
x[1] = 2.2
y[1] (analytic) = 1.6003700000525310174418340844513
y[1] (numeric) = 1.6003700000524449976647129019115
absolute error = 8.60197771211825398e-14
relative error = 5.3749931027424283906659707387082e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 2.83
Order of pole = 0.7924
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
memory used=41.9MB, alloc=4.5MB, time=4.35
x[1] = 2.3
y[1] (analytic) = 1.65703706182307134870637505503
y[1] (numeric) = 1.6570370618229851644481960205557
absolute error = 8.61842581790344743e-14
relative error = 5.2011062495014201100344068338923e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 2.4
y[1] (analytic) = 1.7135309334238865132119933363154
y[1] (numeric) = 1.7135309334238001885160295902036
absolute error = 8.63246959637461118e-14
relative error = 5.0378253628171560721889170876994e-12 %
Correct digits = 13
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
memory used=45.7MB, alloc=4.6MB, time=4.74
x[1] = 2.5
y[1] (analytic) = 1.7698603854199589820629240910936
y[1] (numeric) = 1.769860385419872536970396590062
absolute error = 8.64450925275010316e-14
relative error = 4.8842887969939517058081821229544e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 2.6
y[1] (analytic) = 1.8260335385675518573165010099528
y[1] (numeric) = 1.826033538567465308610076806358
absolute error = 8.65487064242035948e-14
relative error = 4.7397106677513432874171514202247e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 3.121
Order of pole = 0.7898
Complex 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] = 2.7
y[1] (analytic) = 1.8820579264623169828292676921494
y[1] (numeric) = 1.8820579264622303446247184515514
absolute error = 8.66382045492405980e-14
relative error = 4.6033760880088027939706319316866e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
memory used=49.5MB, alloc=4.6MB, time=5.13
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 2.8
y[1] (analytic) = 1.9379405508148655122180207950299
y[1] (numeric) = 1.9379405508147787964377262831156
absolute error = 8.67157802945119143e-14
relative error = 4.4746357290500810742922680330099e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 4.714
Order of pole = 1.101
Complex 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] = 2.9
y[1] (analytic) = 1.9936879303710401970739121219138
y[1] (numeric) = 1.9936879303709534138278530152963
absolute error = 8.67832460591066175e-14
relative error = 4.3529002075543291266983638416307e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 7.65
Order of pole = 2.011
Complex estimate of poles used
memory used=53.4MB, alloc=4.6MB, time=5.53
Complex estimate of poles used
Complex estimate of poles used
x[1] = 3
y[1] (analytic) = 2.0493061443340548456976226184613
y[1] (numeric) = 2.0493061443339680035915155149892
absolute error = 8.68421061071034721e-14
relative error = 4.2376345938944029436402731718058e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 1.188
Order of pole = 0.5348
Complex estimate of poles used
Complex 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] = 3.1
y[1] (analytic) = 2.1048008710137486564058064050106
y[1] (numeric) = 2.1048008710136617627915375466856
absolute error = 8.68936142688583250e-14
relative error = 4.1283532074465172483430661893842e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 0.9619
Order of pole = 0.5219
Complex estimate of poles used
memory used=57.2MB, alloc=4.6MB, time=5.92
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 3.2
y[1] (analytic) = 2.1601774223189895530550478341365
y[1] (numeric) = 2.1601774223189026142351649146254
absolute error = 8.69388198829195111e-14
relative error = 4.0246147832426243361943151748205e-12 %
Correct digits = 13
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
x[1] = 3.3
y[1] (analytic) = 2.2154407746184476386158535973822
y[1] (numeric) = 2.2154407746183606600112918370212
absolute error = 8.69786045617603610e-14
relative error = 3.9260180438242667784508553005934e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=61.0MB, alloc=4.6MB, time=6.32
Complex 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] = 3.4
y[1] (analytic) = 2.270595596419131558902486538415
y[1] (numeric) = 2.2705955964190445451907288606153
absolute error = 8.70137117576777997e-14
relative error = 3.8321976795385209504971347503106e-12 %
Correct digits = 13
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
x[1] = 3.5
y[1] (analytic) = 2.3256462732485114210044978636711
y[1] (numeric) = 2.3256462732484243762338480434299
absolute error = 8.70447706498202412e-14
relative error = 3.7428207226129140678903163510808e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
memory used=64.8MB, alloc=4.6MB, time=6.72
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 3.6
y[1] (analytic) = 2.3805969300725563492610051303924
y[1] (numeric) = 2.3805969300724692769454756898574
absolute error = 8.70723155294405350e-14
relative error = 3.6575832905399372288194826464467e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
x[1] = 3.7
y[1] (analytic) = 2.4354514515368317450717980878139
y[1] (numeric) = 2.4354514515367446482701988303329
absolute error = 8.70968015992574810e-14
relative error = 3.5762076696005246281500201576599e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=68.6MB, alloc=4.6MB, time=7.11
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 3.8
y[1] (analytic) = 2.4902135002795053648859274514507
y[1] (numeric) = 2.4902135002794182462680242044219
absolute error = 8.71186179032470288e-14
relative error = 3.4984397078189762496005421313659e-12 %
Correct digits = 13
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
x[1] = 3.9
y[1] (analytic) = 2.5448865335325435022505521247499
y[1] (numeric) = 2.5448865335324563641526022414873
absolute error = 8.71380979498832626e-14
relative error = 3.4240464870128151021723748514359e-12 %
Correct digits = 13
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
memory used=72.4MB, alloc=4.6MB, time=7.53
Complex estimate of poles used
Complex estimate of poles used
x[1] = 4
y[1] (analytic) = 2.5994738181995926360154858944913
y[1] (numeric) = 2.5994738181995054804870124303551
absolute error = 8.71555284734641362e-14
relative error = 3.3528142450701215667163080863200e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 3.536
Order of pole = 0.7066
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
x[1] = 4.1
y[1] (analytic) = 2.6539784445752622146890724436479
y[1] (numeric) = 2.6539784445751750435323861718407
absolute error = 8.71711566862718072e-14
relative error = 3.2845465216362191424469079732992e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
memory used=76.2MB, alloc=4.6MB, time=7.96
Complex estimate of poles used
x[1] = 4.2
y[1] (analytic) = 2.708403338850112444200878333896
y[1] (numeric) = 2.7084033388500252590045756810967
absolute error = 8.71851963026527993e-14
relative error = 3.2190625026946096747865139691734e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
Complex estimate of poles used
x[1] = 4.3
y[1] (analytic) = 2.7627512745280797404999362051485
y[1] (numeric) = 2.7627512745279925426673762578879
absolute error = 8.71978325599472606e-14
relative error = 3.1561955418821401789627131963580e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=80.1MB, alloc=4.6MB, time=8.39
Complex estimate of poles used
x[1] = 4.4
y[1] (analytic) = 2.8170248828679047690038227601234
y[1] (numeric) = 2.8170248828678175597774057656895
absolute error = 8.72092264169944339e-14
relative error = 3.0957918386652684487645574332065e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 3.026
Order of pole = 0.6315
Complex estimate of poles used
Complex 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] = 4.5
y[1] (analytic) = 2.8712266624470000775574273699597
y[1] (numeric) = 2.8712266624469128580393513727979
absolute error = 8.72195180759971618e-14
relative error = 3.0377092556553655697159086299233e-12 %
Correct digits = 13
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
x[1] = 4.6
memory used=83.9MB, alloc=4.6MB, time=8.82
y[1] (analytic) = 2.9253589879348027185376574425498
y[1] (numeric) = 2.9253589879347154897077116555873
absolute error = 8.72288299457869625e-14
relative error = 2.9818162593223251148990357881065e-12 %
Correct digits = 13
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
x[1] = 4.7
y[1] (analytic) = 2.9794241181527477951864092781936
y[1] (numeric) = 2.9794241181526605579172668591449
absolute error = 8.72372691424190487e-14
relative error = 2.9279909701646110877664715829508e-12 %
Correct digits = 13
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
x[1] = 4.8
y[1] (analytic) = 3.0334242034893580783237709710155
y[1] (numeric) = 3.0334242034892708333941656384962
absolute error = 8.72449296053325193e-14
relative error = 2.8761203100105281325419949936902e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
NO POLE
Complex estimate of poles used
memory used=87.7MB, alloc=4.6MB, time=9.25
Complex estimate of poles used
Complex estimate of poles used
Complex estimate of poles used
Complex 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] = 4.9
y[1] (analytic) = 3.0873612927313928704756583488807
y[1] (numeric) = 3.0873612927313056185817652501906
absolute error = 8.72518938930986901e-14
relative error = 2.8260992355678209940610305570108e-12 %
Correct digits = 13
h = 0.005
TOP MAIN SOLVE Loop
Complex estimate of poles used
Radius of convergence = 2.842
Order of pole = 0.5986
Complex estimate of poles used
Complex estimate of poles used
x[1] = 5
y[1] (analytic) = 3.1412373393653841840133718603913
y[1] (numeric) = 3.1412373393652969257786605370968
absolute error = 8.72582347113232945e-14
relative error = 2.7778300486187345312074868980797e-12 %
Correct digits = 13
h = 0.005
Finished!
diff ( y , x , 1 ) = (0.1 * x + 0.2) / (0.2 * x + 0.3);
Iterations = 980
Total Elapsed Time = 9 Seconds
Elapsed Time(since restart) = 9 Seconds
Time to Timeout = 9 Minutes 50 Seconds
Percent Done = 100.1 %
> quit
memory used=91.2MB, alloc=4.6MB, time=9.65