# BEGIN OUTFILE1 start_glob_y := proc() global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; local r_order,term_no; r_order := 1; while r_order <= 3 do term_no := 1; while term_no <= 3 - r_order + 1 do glob_y_higher[r_order,term_no] := evalf((glob_y_init[term_no + r_order -1])) * glob_h ^ (term_no-1) / ((factorial_1(term_no-1))); term_no := term_no + 1; od; r_order := r_order + 1; od; end; half_jump_glob_y := proc() global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; local r_order,term_no,term_no_adj; r_order := 3; while r_order >= 1 do term_no := glob_max_taylor - 1; while term_no >= 3 - r_order + 1 do term_no_adj := glob_h/term_no; glob_y_higher[r_order,term_no+1] := evalf(glob_y_higher[r_order+1,term_no]*term_no_adj); term_no := term_no - 1; od; r_order := r_order - 1; od; end; jump_glob_y := proc() local r_order,term_no,term_no_adj,calc_term; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; r_order := 3; while r_order >= 1 do term_no := glob_max_taylor - 1; while term_no >= 1 do term_no_adj := evalf(glob_h/term_no); glob_y_higher[r_order,term_no+1] := evalf(glob_y_higher[r_order+1,term_no]*term_no_adj); term_no := term_no - 1; od; calc_term := 1; term_no := glob_max_taylor; glob_y_sav := 0.0; while term_no >= calc_term do glob_y_sav := evalf(copy(glob_y_sav) + glob_y_higher[r_order,term_no]); term_no := term_no - 1; od; glob_y_higher_work[r_order,calc_term] := evalf(glob_y_sav); r_order := r_order - 1; od; r_order := 3; while r_order >= 1 do term_no := r_order; while term_no >= 1 do glob_y_higher[r_order,term_no] := evalf(glob_y_higher_work[r_order,term_no]); term_no := term_no - 1; od; r_order := r_order - 1; od; end; save_data_to_check_if_done := proc() global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; glob_y_higher_work[3 + 1,1] := glob_y_higher[3 + 1,1]; end; check_if_sub_iter_done := proc() global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; glob_sub_iter_done := 1; if abs(glob_y_higher_work[3 + 1,1] - glob_y_higher[3 + 1,1]) > glob_max_trunc_err/100.0 then glob_sub_iter_done := 0; fi; end; get_norms := proc() local i,tt,m; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; set_z(glob_norms); i := 1; while i <= glob_max_taylor do tt := glob_y_higher[1,i]; m := evalf(abs(tt)); if m > glob_norms[i] then glob_norms[i] := m; fi; i := i + 1; od; end; set_z := proc(a) local i; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; i := 1; while i <= glob_max_terms do a[i] := 0.0; i := i + 1; od; end; set_zz := proc(a,mx_ord) local i,j; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; i := 1; while i <= mx_ord + 1 do j := 1; while j <= glob_max_terms do a[i,j] := 0.0; j := j + 1; od; i := i + 1; od; end; set_z2 := proc(a,start) local i; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; i := start + 1; while i <= glob_max_terms do a[i] := 0.0; i := i + 1 od; end; set_zz2 := proc(a,mx_ord,start) local i,j; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; i := 1; while i <= mx_ord + 1 do j := start; while j <= glob_max_terms do a[i,j] := 0.0; j := j + 1; od; i := i + 1; od; end; ats := proc(m,a,b,j) local i,l,ret,ma; ret := 0.0; if (j <= m) then ma := m + 1; i := j; while (i <= m) do l := ma - i; ret := ret + a[i]*b[l]; i := i + 1; od; fi; return(ret); end; att := proc(m,a,b,j) local i,l,ret,ma,al; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; ret := 0.0; if (j <= m) then ma := m + 2; i := j; while (i <= m) do l := ma - i; al := (l - 1); ret := ret + a[i]*b[l]*al; i := i + 1; od; ret := ret / m; fi; return (ret); end; factorial_2 := proc(m,n) global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; if m = 0 and n = 0 then return (1.0); fi; if m = 0 then return (factorial_1(n)); fi; if m = n then return (1.0); fi; return(wrk_prod(n,n-m)); end; const := proc(con) local ret; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; ret := Array(1..(glob_max_terms + 1)); set_z(ret); ret[1] := (con); return(ret); end; atomall := proc() local k; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; k := 1; while k <= 3 do glob_y[k] := glob_y_higher[1,k]; k := k + 1; od; # END OUTFILE1 # BEGIN ATOMHDR1 # emit pre diff $eq_no = 1 i = 1 glob_tmp1[1] := glob_y[1 + 1]/(glob_h ^ 1) * (factorial_2(0,1)); # emit pre mult $eq_no = 1 i = 1 glob_tmp2[1] := (glob_m1[1]*(glob_tmp1[1])); # emit pre add $eq_no = 1 i = 1 glob_tmp3[1] := (const(0.0)[1] + (glob_tmp2[1])); # emit pre assign $eq_no = 1 i = 1 $min_hdrs = 6 # END ATOMHDR1 # BEGIN ATOMHDR2 # emit pre diff $eq_no = 1 i = 2 glob_tmp1[2] := glob_y[1 + 2]/(glob_h ^ 1) * factorial_2(2 - 1 ,1 + 2 - 1); # emit pre mult $eq_no = 1 i = 2 glob_tmp2[2] := ats(2,glob_m1,glob_tmp1,1); # emit pre add $eq_no = 1 i = 2 glob_tmp3[2] := (const(0.0)[2] + (glob_tmp2[2])); # END ATOMHDR2 # BEGIN ATOMHDR3 # emit pre diff $eq_no = 1 i = 3 glob_tmp1[3] := glob_y[1 + 3]/(glob_h ^ 1) * factorial_2(3 - 1 ,1 + 3 - 1); # emit pre mult $eq_no = 1 i = 3 glob_tmp2[3] := ats(3,glob_m1,glob_tmp1,1); # emit pre add $eq_no = 1 i = 3 glob_tmp3[3] := (const(0.0)[3] + (glob_tmp2[3])); # END ATOMHDR3 # BEGIN ATOMHDR4 # emit pre diff $eq_no = 1 i = 4 glob_tmp1[4] := glob_y[1 + 4]/(glob_h ^ 1) * factorial_2(4 - 1 ,1 + 4 - 1); # emit pre mult $eq_no = 1 i = 4 glob_tmp2[4] := ats(4,glob_m1,glob_tmp1,1); # emit pre add $eq_no = 1 i = 4 glob_tmp3[4] := (const(0.0)[4] + (glob_tmp2[4])); if 7 <= glob_max_taylor then glob_y[7] := glob_tmp3[4] * (glob_h ^ 3) / factorial_2(3,6); fi; # END ATOMHDR4 # BEGIN ATOMHDR5 # emit pre diff $eq_no = 1 i = 5 glob_tmp1[5] := glob_y[1 + 5]/(glob_h ^ 1) * factorial_2(5 - 1 ,1 + 5 - 1); # emit pre mult $eq_no = 1 i = 5 glob_tmp2[5] := ats(5,glob_m1,glob_tmp1,1); # emit pre add $eq_no = 1 i = 5 glob_tmp3[5] := (const(0.0)[5] + (glob_tmp2[5])); if 8 <= glob_max_taylor then glob_y[8] := glob_tmp3[5] * (glob_h ^ 3) / factorial_2(4,7); fi; # END ATOMHDR5 # BEGIN ATOMHDR6 # emit pre diff $eq_no = 1 i = 6 glob_tmp1[6] := glob_y[1 + 6]/(glob_h ^ 1) * factorial_2(6 - 1 ,1 + 6 - 1); # emit pre mult $eq_no = 1 i = 6 glob_tmp2[6] := ats(6,glob_m1,glob_tmp1,1); # emit pre add $eq_no = 1 i = 6 glob_tmp3[6] := (const(0.0)[6] + (glob_tmp2[6])); if 9 <= glob_max_taylor then glob_y[9] := glob_tmp3[6] * (glob_h ^ 3) / factorial_2(5,8); fi; # END ATOMHDR6 # BEGIN OUTFILE3 k := 6 + 0 + 1; while k <= glob_max_taylor do # END OUTFILE3 # BEGIN OUTFILE4 # emit diff $eq_no = 1 if k + 1 <= glob_max_taylor then glob_tmp1[k] := glob_y[k + 1] / (glob_h ^ 1) * factorial_2(k - 1,k + 1-1); fi; # emit mult $eq_no = 1 glob_tmp2[k] := ats(k,glob_m1,glob_tmp1,1); # emit add $eq_no = 1 glob_tmp3[k] := (const(0.0)[k]+(glob_tmp2[k])); # emit assign $eq_no = 1 if k > 3 then if k + 3 <= glob_max_taylor then glob_y[k + 3 ] := glob_tmp3[k] * (glob_h ^ 3) / factorial_2( k-1 ,3+k-1); fi; fi; k := k + 1; od; k := 1; while k <= glob_max_taylor do glob_y_higher[3 + 1,k] := glob_tmp3[k]; k := k + 1; od; # END OUTFILE4 # BEGIN OUTFILE5 end; # BEGIN ATS LIBRARY BLOCK ################################################################### ### # # File: Ats.mxt # # Subject: Maple library for a Maple program ("diffeq.mxt") # generated by Sode.rb # to solve a system of ordinary differential # equations with long Taylor series. # # Author: Dennis J. Darland # # Date: June 20, 2008 # ############################################################################ # # Copyright (C) 2008 Dennis J. Darland# # # This program is free software you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software 3) and ((abs(y_higher[1,m]) < glob_small_float) or (abs(y_higher[1,m-1]) < glob_small_float) or (abs(y_higher[1,m-2]) < glob_small_float )) do # if DJD_DEBUG = 1 then # print("radius_real m = " , m); # fi; m := m - 1; od; if (m > 3) then rm0 := y_higher[1,m]/y_higher[1,m-1]; rm1 := y_higher[1,m-1]/y_higher[1,m-2]; hdrc := (m-1)*rm0-(m-2)*rm1; if DJD_DEBUG = 1 then print("re hdrc = " , hdrc); fi; if abs(hdrc) > glob_small_float then rcs := glob_h/hdrc; ord_no := (m-1)*rm0/hdrc -m + 2; ans[1] := rcs; ans[2] := ord_no; if DJD_DEBUG = 1 then print("ans[1] = " , ans[1]); print("ans[2] = " , ans[2]); fi; return(ans); else ans[1] := glob_large_float; ans[2] := glob_large_float; return(ans); fi; else ans[1] := glob_large_float; ans[2] := glob_large_float; return(ans); fi; end; radius_complex := proc(y_higher,order_diff) local ans,n,cnt,m,rm0,rm1,rm2,rm3,rm4,nr1,nr2,dr1,dr2,ds1,ds2,rcs,ord_no,rad_c; global glob_h,glob_large_float, glob_max_taylor, glob_small_float; ## # 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 ans := Array(1..3); n := glob_max_taylor; cnt := 0; while (cnt < 5 and n >= 1) do if ((abs(y_higher[1,n])) > (glob_small_float)) then cnt := cnt + 1; else cnt := 0; fi; n := n - 1; od; m := n + cnt; if (m < 6) then ans[1] := (glob_large_float); ans[2] := (glob_large_float); return(ans); fi; if abs(y_higher[1,m]) >= (glob_large_float) or abs(y_higher[1,m-1]) >=(glob_large_float) or abs(y_higher[1,m-2]) >= (glob_large_float) or abs(y_higher[1,m-3]) >= (glob_large_float) or abs(y_higher[1,m-4]) >= (glob_large_float) or abs(y_higher[1,m-5]) >= (glob_large_float) then ans[1] := (glob_large_float); ans[2] := (glob_large_float); return(ans); else rm0 := (y_higher[1,m])/(y_higher[1,m-1]); rm1 := (y_higher[1,m-1])/(y_higher[1,m-2]); rm2 := (y_higher[1,m-2])/(y_higher[1,m-3]); rm3 := (y_higher[1,m-3])/(y_higher[1,m-4]); rm4 := (y_higher[1,m-4])/(y_higher[1,m-5]); nr1 := (m-1)*rm0 - 2.0*(m-2)*rm1 + (m-3)*rm2; nr2 := (m-2)*rm1 - 2.0*(m-3)*rm2 + (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 abs(nr1 * dr2 - nr2 * dr1) <= glob_small_float or abs(dr1) <= glob_small_float then ans[1] := (glob_large_float); ans[2] := (glob_large_float); return(ans); else #term + dr1*dr2 added to Manuel Prieto's analysis if abs(nr1*dr2 - nr2 * dr1) > glob_small_float then rcs := ((ds1*dr2 - ds2*dr1 +dr1*dr2)/(nr1*dr2 - nr2 * dr1)); # rcs := (ds1*dr2 - ds2*dr1)/(nr1*dr2 - nr2 * dr1) ord_no := (rcs*nr1 - ds1)/(2.0*dr1) -(m)/2.0; if abs(rcs) > glob_small_float then if rcs > 0.0 then rad_c :=sqrt(rcs) * glob_h; else rad_c :=glob_large_float; fi; else rad_c := glob_large_float; ord_no := glob_large_float; fi; else rad_c := glob_large_float; ord_no := glob_large_float; fi; ans[1] := rad_c; ans[2] := ord_no; return(ans); fi; fi; end; which_radius := proc(prev_pole,real_pole,complex_pole,print_flag) local pole; global glob_correct_start_flag, glob_large_float,glob_out_fd; pole := Array(1..3); if (real_pole[1] = (glob_large_float) or real_pole[2] = (glob_large_float)) and (complex_pole[1] <> (glob_large_float)) and (complex_pole[2] <> (glob_large_float)) and (complex_pole[1] > 0.0) and (complex_pole[2] > (0.0)) then pole[1] := complex_pole[1]; pole[2] := complex_pole[2]; if glob_correct_start_flag <> 1 and print_flag = 1 then print("complex estimate of poles used 1"); fprintf(glob_out_fd,"complex estimate of poles used 1\n"); fi; print_pole(pole); return(pole); fi; if (real_pole[1] <> glob_large_float) and (real_pole[2] <> glob_large_float) and (real_pole[1] > 0.0) and (real_pole[2] > 0.0) and ((complex_pole[1] = glob_large_float) or (complex_pole[2] = glob_large_float) or (complex_pole[1] <= 0.0 ) or (complex_pole[2] <= 0.0)) then pole[1] := real_pole[1]; pole[2] := real_pole[2]; if glob_correct_start_flag <> 1 and print_flag = 1 then print("real estimate of pole used 1"); fprintf(glob_out_fd,"real estimate of pole used 1\n"); fi; print_pole(pole); return(pole); fi; if ((real_pole[1] = glob_large_float) or (real_pole[2] = glob_large_float)) and ((complex_pole[1] = glob_large_float) or (complex_pole[2] = glob_large_float)) then pole[1] := glob_large_float; pole[2] := glob_large_float; if glob_correct_start_flag <> 1 and print_flag = 1 then print("NO POLE 1"); fprintf(glob_out_fd,"NO POLE 1\n"); fi; print_pole(pole); return(pole); fi; if (real_pole[1] < complex_pole[1]) and (real_pole[1] > 0.0) and (real_pole[2] > 0.0) then pole[1] := real_pole[1]; pole[2] := real_pole[2]; if glob_correct_start_flag <> 1 and print_flag = 1 then print("real estimate of pole used 2"); fprintf(glob_out_fd,"real estimate of pole used 2\n"); fi; print_pole(pole); return(pole); fi; if (complex_pole[1] <> glob_large_float) and (complex_pole[2] <> glob_large_float) and (complex_pole[1] > 0.0) and (complex_pole[2] > 0.0) then pole[1] := complex_pole[1]; pole[2] := complex_pole[2]; if glob_correct_start_flag <> 1 and print_flag = 1 then print("complex estimate of poles used 2"); fprintf(glob_out_fd,"complex estimate of poles used 2\n"); fi; print_pole(pole); return(pole); fi; pole[1] := glob_large_float; pole[2] := glob_large_float; if glob_correct_start_flag <> 1 and print_flag = 1 then print("NO POLE 2"); fprintf(glob_out_fd,"NO POLE 2\n"); fi; return(pole); end; print_pole := proc(pole) global glob_correct_start_flag, glob_large_float, print_flag,glob_out_fd; if (pole[1] <> glob_large_float) and (glob_correct_start_flag <> 1 and print_flag = 1) then print("radius of convergence = " , pole[1]); print("r_order of singularity = " , pole[2]); fprintf(glob_out_fd,"radius of convergence = %20.16g\n" , pole[1]); fprintf(glob_out_fd,"r_order of singularity = %20.16g\n" , pole[2]); fi; end; compute_order := proc(log10norm, log10abserr, log10relerr) global glob_max_taylor; #* # * this is to determine the 'optimal' degree of the taylor expansion. # * # * parameters: # * log10norm: base 10 log of the norm of the initial condition # * log10abserr: base-10 log of the absolute error required # * log10relerr: base-10 log of the relative error required # * # * returned value: 'optimal' degree. # # When I studied it I couldn't really understand Angel Jorba and Maorong Zou 's # method for determining the r_order. I don't adjust it anymore. # It will be set to the problem file max_terms which defaults to 30 which seems # a good value most of the time. # I've retained glob_max_taylor in case I find a better way. return(glob_max_taylor); end; compute_step_size := proc(x,print_flag,log10norm,log10abserr,log10relerr) # # * it looks for a step size to keep the truncation error (both # * absolute and relative) below the absolute error tol, and also tries # * to reduce cancellations of big terms in the summation of the taylor # * series. # * # * WARNING: there are two versions of this function. this is the one # * that uses MY_FLOAT variables for all the computations # * local dh,h,hs,hd,z,hdj,flag,j,r,a,hd2; global glob_adjusted_h_last_time,glob_correct_start_flag, glob_h, glob_hmax, glob_hmin_init,glob_max_taylor, glob_norms, glob_reached_optimal_h, DJD_DEBUG; #local dh,hd,j,aaa,r2,cc,hd2,i #local h,r,z,hs cs_info("begin compute step size"); if DJD_DEBUG = 1 then print("compute_step_size last time := " , glob_adjusted_h_last_time , "optimal := " , glob_reached_optimal_h , "start_flag := " , glob_correct_start_flag); fi; if glob_adjusted_h_last_time = 1 or glob_reached_optimal_h = 1 then glob_adjusted_h_last_time := 0; return glob_h; fi; # we will use the array norms to store the norms of the normalized # derivatives #*/ #/* # we compute the step size according to the first algorithm #*/ cs_info("begin compute step size A"); get_norms(); dh:=compute_step_size_A(x,print_flag,log10norm,log10abserr,log10relerr); cs_info("After compute step size A"); h := (dh); hs := (h); #/* # next loop runs over the norms of the taylor coefficients, checking # whether the summation of the taylor series (from low to high r_order) # involves increasing terms. in this case, it reduces the stepsize # accordingly. #*/ hd := h; z := glob_norms[1]*hd; hdj := hd; flag := 0; j := 1; while j<= glob_max_taylor-1 do hdj := hdj * hd; r:= glob_norms[j+1]*hdj; if (r <= z) then j := j + 1; next; fi; # reduce hd (& hdj) flag := 1; hdj := z/glob_norms[j+1]; a := (r/z) ^ j; hd2 := hd; hd:=hd2/a; # z := (r) j := j + 1; od; if hd < (glob_hmin_init) then return (glob_hmin_init); fi; if hd > (glob_hmax) then return (glob_hmax); fi; return (hd); end; estimate_trunc_err := proc() local d1,d2,d3,d4; global glob_max_order, glob_max_taylor, glob_norms; get_norms(); d1 := glob_norms[glob_max_taylor-glob_max_order]; d2 := glob_norms[glob_max_taylor-1-glob_max_order]; d3 := glob_norms[glob_max_taylor-2-glob_max_order]; d4 := glob_norms[glob_max_taylor-3-glob_max_order]; return (abs(d1)+abs(d2)+abs(d3)+abs(d4)); end; compute_step_size_A := proc(x,print_flag,log10norm,log10abserr,log10relerr) local d1,d2,d3,d4,log10d1,log10d2,log10d3,log10d4,r,log10ro4,log10ro3, log10ro2,log10ro1,a_log10ro,size_it; global glob_unchanged_h_cnt, glob_reached_optimal_h, glob_optimal_clock_start_sec, glob_optimal_start, glob_h, glob_hmax, glob_hmin_init, glob_log10abserr, glob_max_order, glob_max_taylor, glob_no_eqs, glob_norms, glob_small_float; if log10abserr < (log10relerr + log10norm) then log10abserr := log10relerr + log10norm; fi; d1 := max(glob_norms[glob_max_taylor-glob_max_order],glob_small_float); d2 := max(glob_norms[glob_max_taylor-1-glob_max_order],glob_small_float); d3 := max(glob_norms[glob_max_taylor-2-glob_max_order],glob_small_float); d4 := max(glob_norms[glob_max_taylor-3-glob_max_order],glob_small_float); log10d1 := evalf(log10(d1)); log10d2 := evalf(log10(d2)); log10d3 := evalf(log10(d3)); log10d4 := evalf(log10(d4)); r := log10abserr; log10ro4 := evalf(r - log10d4)/evalf(glob_max_taylor-4-glob_max_order); log10ro3 := evalf(r - log10d3)/evalf(glob_max_taylor-3-glob_max_order); log10ro2 := evalf(r - log10d2)/evalf(glob_max_taylor-2-glob_max_order); log10ro1 := evalf(r - log10d1)/evalf(glob_max_taylor-1-glob_max_order); a_log10ro := max(log10ro1,log10ro2,log10ro3,log10ro4); size_it := 10.0 ^ a_1og10ro; cs_info("After size_it computed"); if glob_hmax < size_it and glob_h <> glob_hmin_init then glob_unchanged_h_cnt := 0; glob_reached_optimal_h := 1; glob_optimal_clock_start_sec := time(); glob_optimal_start := x[1]; if glob_hmax < size_it then size_it:= glob_hmax; fi; fi; if glob_h * 2.0 < size_it then print("glob_h suggested at " , (glob_h * 2.0)); print("Requested accuracy was glob_log10abserr = " , glob_log10abserr); cs_info("glob_h suggested "); return (2.0 * glob_h); fi; return(size_it); end; adjust_step := proc(x,y_higher,hold,hnew,order_diff) local r_order,term_no,accum,r; global glob_hmin_init, glob_log10abserr, glob_max_terms; if hnew < glob_hmin_init then hnew := glob_hmin_init; fi; if hold >= hnew * 2.0 then print("Adjusting glob_h from " , hold , " to " , hnew); print("Requested accuracy was glob_log10abserr = " , glob_log10abserr); #glob_x[2]:= hnew; accum:= 1.0; r := hnew/hold; term_no := 1; while term_no <= glob_max_terms do r_order := 1; while r_order <= order_diff + 1 do y_higher[r_order,term_no] := y_higher[r_order,term_no] * accum; r_order := r_order + 1; od; accum := accum*r; term_no := term_no + 1; od; return (hnew); else return (hold); fi; end; wrk_prod := proc(f,n) local p,i,r; p := f; r := f - 1; i := n; while (i >= 1 and r >= 1) do p := p * r; r := r - 1; i := i - 1; od; return (p); end; chk_data := proc() local errflag; global glob_max_terms, glob_max_taylor, glob_look_poles, glob_adjust_h, glob_max_iter; #local errflag errflag := 0; if glob_max_terms < 15 or glob_max_terms > 512 then print("Illegal max_terms: " , glob_max_terms , " -- Using 30"); glob_max_terms := 30; glob_max_taylor := 30; fi; if (glob_look_poles <> 0 and glob_look_poles <> 1) then print("Illegal look_poles"); errflag := 1; fi; if (glob_adjust_h <> 0 and glob_adjust_h <> 1) then print("Illegal adjust_h "); errflag := 1; fi; if (glob_look_poles = 1 and glob_adjust_h = 0) then print("Illegal look_poles or adjust_h"); errflag := 1; fi; if (glob_max_iter < 10) then print("Illegal max_iter"); errflag := 1; fi; if errflag = 1 then quit; fi; end; comp_expect_sec := proc(t_end2,t_start2,t2,clock_sec) local ms2,r,sub1,sub2,sec_left; global glob_small_float; ms2 := clock_sec; sub1 := (t_end2-t_start2); sub2 := (t2-t_start2); if sub1 = 0.0 then sec_left := 0.0 else if abs(sub2) > glob_small_float then r := (sub1/sub2); sec_left := r * ms2 - ms2; else sec_left := -1.0; fi; fi; return sec_left; end; comp_percent := proc(t_end2,t_start2,t2) local r,sub1,sub2; global glob_small_float; sub1 := (t_end2-t_start2); sub2 := (t2-t_start2); if abs(sub2) > glob_small_float then r := 100.0*(sub2/sub1); else r := 0.0; fi; return round(r); end; time_elapsed := proc(clock_sec_v) local sec,hrs; sec := 1.0*clock_sec_v; hrs := sec/(3600.0); return hrs; end; sec_to_display := proc(sec) global CENTURIES_IN_MILLINIUM, DAYS_IN_YEAR, HOURS_IN_DAY, MIN_IN_HOUR, SEC_IN_MIN, YEARS_IN_CENTURY; local sec_in_millinium,milliniums,millinium_int,centuries,cent_int, years,years_int,days,days_int,hours,hours_int,minutes,minutes_int,seconds, sec_int,str; if sec < 0.0 then return " Unknown" fi; sec_in_millinium := SEC_IN_MIN * MIN_IN_HOUR * HOURS_IN_DAY * DAYS_IN_YEAR * YEARS_IN_CENTURY * CENTURIES_IN_MILLINIUM; milliniums := sec/sec_in_millinium; millinium_int := trunc(milliniums); centuries := (milliniums - millinium_int)* CENTURIES_IN_MILLINIUM; cent_int := trunc(centuries); years := (centuries-cent_int) * YEARS_IN_CENTURY; years_int := trunc(years); days := (years - years_int) * DAYS_IN_YEAR; days_int := trunc(days); hours := (days - days_int) * HOURS_IN_DAY; hours_int := trunc(hours); minutes := (hours - hours_int) * MIN_IN_HOUR; minutes_int := trunc(minutes); seconds := (minutes - minutes_int) * SEC_IN_MIN; sec_int := trunc(seconds); str := " "; if millinium_int > 0 then str := " " || str || millinium_int || " Millinia "; fi; if cent_int > 0 then str := " " || str || cent_int || " Centuries "; fi; if years_int > 0 then str := " " || str || years_int || " Years "; fi; if days_int > 0 then str := " " || str || days_int || " Days "; fi; if hours_int > 0 then str := " " || str || hours_int || " Hours "; fi; if minutes_int > 0 then str := " " || str || minutes_int || " Minutes "; fi; if sec_int > 0 then str := " " || str || sec_int || " Seconds "; fi; return str; end; cs_info := proc(str) global DJD_DEBUG, glob_correct_start_flag, glob_h, glob_reached_optimal_h; if DJD_DEBUG = 1 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; factorial_1 := proc(n) return (n!); end; # END ATS LIBARY BLOCK # BEGIN USER DEF BLOCK exact_soln_glob_y := proc(x) return(sin(x)); end; exact_soln_glob_yp := proc(x) return(cos(x)); end; exact_soln_glob_ypp := proc(x) return(-sin(x)); end; # END USER DEF BLOCK # END OUTFILE5 # BEGIN OUTFIEMAIN local percent_done; global glob_y, glob_tmp0, glob_tmp1, glob_tmp2, glob_m1, glob_tmp3, glob_x, glob_y_init, glob_y_higher, glob_y_higher_work, glob_y_sav, DJD_DEBUG, DJD_DEBUG2, SEC_IN_MIN, MIN_IN_HOUR, HOURS_IN_DAY, DAYS_IN_YEAR, YEARS_IN_CENTURY, CENTURIES_IN_MILLINIUM, glob_out_fd, glob_progress_fd, glob_adjust_h, glob_adjusted_h_last_time, glob_almost_1, glob_clock_sec, glob_clock_start_sec, glob_correct_start_flag, glob_disp_incr, glob_h, glob_hmax, glob_hmin, glob_hmin_init, glob_large_float, glob_last_good_h, glob_look_poles, glob_log10_abserr, glob_log10_relerr, glob_max_hours, glob_max_iter, glob_max_order, glob_max_rel_trunc_err, glob_max_taylor, glob_max_terms, glob_max_trunc_err, glob_no_eqs, glob_norms, glob_old_log10eps, glob_old_nt, glob_optimal_clock_start_sec, glob_optimal_start, glob_progress_ind, glob_small_float, glob_smallish_float, glob_tot_order, glob_unchanged_h_cnt, glob_sub_iter, glob_sub_iter_done, glob_warned, glob_warned2; glob_orig_start_sec := time(); # printlevel := 1; glob_out_fd := fopen(`c:\\cygwin\\home\\dennis\\src\\mine\\MapleSode\\TXT\\diff2.ode.out.txt`,WRITE,TEXT); glob_progress_fd := fopen(`c:\\cygwin\\home\\dennis\\src\\mine\\MapleSode\\progress.txt`,WRITE,TEXT); DJD_DEBUG := 1; DJD_DEBUG2 := 1; SEC_IN_MIN := 60.0; MIN_IN_HOUR := 60.0; HOURS_IN_DAY := 24.0; DAYS_IN_YEAR := 365.0; YEARS_IN_CENTURY := 100.0; CENTURIES_IN_MILLINIUM := 10.0; # trace(start_glob_y); # trace(atomall); # trace(half_jump_glob_y); # trace(jump_glob_y); MAX_UNCHANGED := 0; glob_curr_iter_when_opt := 0; glob_max_terms := 30; glob_max_taylor := 30; print_flag := 1; glob_min_terms := 30; glob_max_order := 3; glob_no_eqs = 1; glob_tot_order = 4; # BEGIN CONST BLOCK glob_max_terms := 30; # END CONST BLOCK iter := -1; glob_norms := Array(1..(glob_max_terms + 1)); glob_max_iter := 50000; glob_max_hours := 0.0; glob_max_minutes := 15.0; glob_y := Array(1..(glob_max_terms + 1)); set_z(glob_y); glob_tmp0 := Array(1..(glob_max_terms + 1)); set_z(glob_tmp0); glob_tmp1 := Array(1..(glob_max_terms + 1)); set_z(glob_tmp1); glob_tmp2 := Array(1..(glob_max_terms + 1)); set_z(glob_tmp2); glob_m1 := Array(1..(glob_max_terms + 1)); set_z(glob_m1); glob_tmp3 := Array(1..(glob_max_terms + 1)); set_z(glob_tmp3); glob_x := Array(1..(glob_max_terms + 1)); set_z(glob_x); glob_y_init := Array(1..(glob_max_terms + 1)); glob_y_higher := Array(1..4,1..(glob_max_terms+1)); glob_y_higher_work := Array(1..4,1..(glob_max_terms+1)); print("##############ECHO OF PROBLEM#################"); fprintf(glob_out_fd,"##############ECHO OF PROBLEM#################\n"); print("##############diff2.ode#################"); fprintf(glob_out_fd,"##############diff2.ode#################\n"); print("diff ( y , x , 3 ) = m1 * diff ( y , x , 1 ) ;"); fprintf(glob_out_fd,"diff ( y , x , 3 ) = m1 * diff ( y , x , 1 ) ;\n"); print("!"); fprintf(glob_out_fd,"!\n"); print("glob_max_terms := 30;"); fprintf(glob_out_fd,"glob_max_terms := 30;\n"); print("!"); fprintf(glob_out_fd,"!\n"); print("glob_x_start := 0.1;"); fprintf(glob_out_fd,"glob_x_start := 0.1;\n"); print("glob_x_end := 0.6 ;"); fprintf(glob_out_fd,"glob_x_end := 0.6 ;\n"); print("glob_y_init[1] := exact_soln_glob_y(glob_x_start);"); fprintf(glob_out_fd,"glob_y_init[1] := exact_soln_glob_y(glob_x_start);\n"); print("glob_y_init[2] := exact_soln_glob_yp(glob_x_start);"); fprintf(glob_out_fd,"glob_y_init[2] := exact_soln_glob_yp(glob_x_start);\n"); print("glob_y_init[3] := exact_soln_glob_ypp(glob_x_start);"); fprintf(glob_out_fd,"glob_y_init[3] := exact_soln_glob_ypp(glob_x_start);\n"); print("glob_hmax := 0.0001;"); fprintf(glob_out_fd,"glob_hmax := 0.0001;\n"); print("glob_adjust_h := 1;"); fprintf(glob_out_fd,"glob_adjust_h := 1;\n"); print("glob_look_poles := 0;"); fprintf(glob_out_fd,"glob_look_poles := 0;\n"); print("glob_max_iter := 400;"); fprintf(glob_out_fd,"glob_max_iter := 400;\n"); print("glob_log10relerr := -8.0;"); fprintf(glob_out_fd,"glob_log10relerr := -8.0;\n"); print("glob_log10abserr := -8.0;"); fprintf(glob_out_fd,"glob_log10abserr := -8.0;\n"); print("!"); fprintf(glob_out_fd,"!\n"); print("exact_soln_glob_y := proc(x)"); fprintf(glob_out_fd,"exact_soln_glob_y := proc(x)\n"); print(" return(sin(x));"); fprintf(glob_out_fd," return(sin(x));\n"); print("end;"); fprintf(glob_out_fd,"end;\n"); print("exact_soln_glob_yp := proc(x)"); fprintf(glob_out_fd,"exact_soln_glob_yp := proc(x)\n"); print(" return(cos(x));"); fprintf(glob_out_fd," return(cos(x));\n"); print("end;"); fprintf(glob_out_fd,"end;\n"); print("exact_soln_glob_ypp := proc(x)"); fprintf(glob_out_fd,"exact_soln_glob_ypp := proc(x)\n"); print(" return(-sin(x));"); fprintf(glob_out_fd," return(-sin(x));\n"); print("end;"); fprintf(glob_out_fd,"end;\n"); print(""); fprintf(glob_out_fd,"\n"); print("!"); fprintf(glob_out_fd,"!\n"); print("#######END OF ECHO OF PROBLEM#################"); fprintf(glob_out_fd,"#######END OF ECHO OF PROBLEM#################\n"); set_z(glob_norms); set_z(glob_y_init); set_zz(glob_y_higher,3); set_zz(glob_y_higher_work,3); glob_correct_start_flag := 1; glob_unchanged_h_cnt := 0; glob_adjust_h := 1; Digits := 20; glob_look_poles := 0; glob_max_sub_iter := 256; glob_min_sub_iter := 2; glob_warned := 0; glob_warned2 := 0; glob_small_float := 1.0e-200; glob_smallish_float := 1.0e-10; glob_large_float := 1.0e100; glob_almost_1 := 0.99; glob_progress_ind := 1; glob_log10abserr := -8; glob_log10relerr := -8; glob_hmax := 0.01; glob_old_log10eps := -1.0; glob_old_nt := 0; glob_adjusted_h_last_time := 0; glob_reached_optimal_h := 0; # djd_cnt := 0; # BEGIN INPUT BLOCK glob_x_start := 0.1; glob_x_end := 0.6 ; glob_y_init[1] := exact_soln_glob_y(glob_x_start); glob_y_init[2] := exact_soln_glob_yp(glob_x_start); glob_y_init[3] := exact_soln_glob_ypp(glob_x_start); glob_hmax := 0.0001; glob_adjust_h := 1; glob_look_poles := 0; glob_max_iter := 400; glob_log10relerr := -8.0; glob_log10abserr := -8.0; # END INPUT BLOCK glob_hmin_init := evalf(10.0 ^ (glob_log10abserr/4.0)); glob_max_trunc_err := evalf(10.0 ^ glob_log10abserr); glob_max_rel_trunc_err := evalf(10.0 ^ glob_log10relerr); hnew := glob_hmin_init; glob_last_good_h := hnew; glob_adjust_h := 1; glob_progress_ind := 1; glob_max_sec := evalf(60.0 * (glob_max_minutes + 60.0 * glob_max_hours)); cs_info("before cs loop"); while glob_correct_start_flag <> 0 do # print("DENNIS GOT HERE!!!"); cs_info("start loop"); glob_h := evalf(hnew); glob_x[1] := glob_x_start; glob_x[2] := (glob_h); glob_hsav := evalf(glob_h); glob_max_taylor := glob_max_terms; glob_max_taylor_SAV := glob_min_terms; chk_data(); log10normmin := evalf(log10(glob_large_float)); start_glob_y(); if abs(glob_y_higher[1,1]) > glob_small_float then tmp := abs(glob_y_higher[1,1]); log10norm := evalf(log10(tmp)); if log10norm < log10normmin then log10normmin := log10norm; fi; fi; if glob_h = glob_hmin_init and iter = -1 then print("Testing for glob_h"); print("Finding the largest glob_h to give the accuracy needed saves time!"); print("You must restart if you increase glob_h as recalculating the Taylor series terms results in HUGE errors!"); print("We start with an glob_h based on the requested error glob_h :=" , glob_h); fi; if glob_h = glob_hmin_init and iter = -1 then print_alot(1,glob_x,glob_y_higher,3,"glob_y","glob_x") fi; glob_clock_start_sec := time(); cs_info("main loop loop"); glob_clock_sec := time(); current_iter := 0; while (glob_x[1] <= (glob_x_end ) and iter < glob_max_iter and (glob_clock_sec - glob_orig_start_sec) < glob_max_sec) do iter := iter + 1; current_iter := current_iter + 1; atomall(); half_jump_glob_y(); cs_info("before progress report"); if ((iter mod glob_progress_ind) = 0) then clock_sec1 := time(); total_clock_sec := clock_sec1 - glob_orig_start_sec; total_disp := sec_to_display(total_clock_sec); clock_sec := clock_sec1 - glob_clock_start_sec; elapsed_disp := sec_to_display(clock_sec); left_sec := glob_max_sec + glob_orig_start_sec - clock_sec1; time_till_disp := sec_to_display(left_sec); if glob_reached_optimal_h = 1 then expect_sec := comp_expect_sec(glob_x_end,glob_x_start,glob_x[1] , clock_sec); opt_clock_sec := clock_sec1 - glob_optimal_clock_start_sec; optimal_expect_sec := comp_expect_sec(glob_x_end,glob_optimal_start,glob_x[1] , opt_clock_sec); percent_done := comp_percent(glob_x_end,glob_optimal_start,glob_x[1]); optimal_expect_disp := sec_to_display(optimal_expect_sec); else expect_sec := comp_expect_sec(glob_x_end,glob_x_start,glob_x[1] + glob_h , clock_sec); optimal_expect_disp := " Unknown"; fi;; expect_disp := sec_to_display(expect_sec); print(" " || "diff2.ode" || " glob_h = " || glob_h || " glob_max_taylor = " || glob_max_taylor || " Iterations = " || iter || " Current Iterations = " || current_iter); print(" " || "Total Elapsed Time = " || total_disp); print(" " || "Elapsed Time(since restart)= " || elapsed_disp); print(" " || "Expected Time Remaining = " || expect_disp); print(" " || "Optimized Time Remaining = " || optimal_expect_disp); print(" " || "Time to Timeout = " || time_till_disp); print(" " || "Percent Done = " || percent_done); fprintf(glob_out_fd," " || "diff2.ode glob_h = %g glob_max_taylor = %d Iterations = %d Current Iterations = %d\n",glob_h,glob_max_taylor,iter,current_iter); fprintf(glob_out_fd," " || "Total Elapsed Time = " || total_disp || "\n"); fprintf(glob_out_fd," " || "Elapsed Time(since restart)= " || elapsed_disp || "\n"); fprintf(glob_out_fd," " || "Expected Time Remaining = " || expect_disp || "\n"); fprintf(glob_out_fd," " || "Optimized Time Remaining = " || optimal_expect_disp || "\n"); fprintf(glob_out_fd," " || "Time to Timeout = " || time_till_disp || "\n"); fprintf(glob_progress_fd," " || "diff2.ode glob_h = %g glob_max_taylor = %d Iterations = %d Current Iterations = %d\n",glob_h,glob_max_taylor,iter,current_iter); fprintf(glob_progress_fd," " || "Total Elapsed Time = " || total_disp || "\n"); fprintf(glob_progress_fd," " || "Elapsed Time(since restart)= " || elapsed_disp || "\n"); fprintf(glob_progress_fd," " || "Expected Time Remaining = " || expect_disp || "\n"); fprintf(glob_progress_fd," " || "Optimized Time Remaining = " || optimal_expect_disp || "\n"); fprintf(glob_progress_fd," " || "Time to Timeout = " || time_till_disp || "\n"); fi; fflush(glob_progress_fd); fflush(glob_out_fd); pole := Array(1..3); pole[1] := glob_large_float; pole[2] := glob_xlarge_float; if glob_look_poles = 1 then real_pole := radius_real(glob_y_higher,3); complex_poles := radius_complex(glob_y_higher,3); pole := which_radius(pole,real_pole,complex_poles,print_flag); fi; if glob_reached_optimal_h = 0 then normmax := glob_small_float; if abs(glob_y_higher[1,1]) > glob_small_float then tmp := abs(glob_y_higher[1,1]); if tmp < normmax then normmax := tmp; fi; fi; fi; if glob_look_poles = 1 and abs(pole[1]) > glob_small_float and (pole[1] <> glob_large_float) then sz2 := pole[1]/10.0; if sz2 < hnew then print("glob_h adjusted to " || sz2 || " due to singularity."); print("Reached Optimal"); if glob_reached_optimal_h = 0 then glob_reached_optimal_h := 1; glob_curr_iter_when_opt := current_iter; glob_correct_start_flag := 0; glob_optimal_clock_start_sec := time(); glob_optimal_start := glob_x[1]; fi; hnew := sz2; fi; fi; if glob_reached_optimal_h = 0 then if glob_max_trunc_err < glob_max_rel_trunc_err * normmax then max_trunc_err := glob_max_rel_trunc_err * normmax / 100.0; else max_trunc_err := glob_max_rel_trunc_err / 100.0; fi; cs_info("BEFORE PROBLEM hnew = " || hnew || " glob_last_good_h = " || glob_last_good_h || " glob_unchanged_h_cnt = " || glob_unchanged_h_cnt || " glob_reached_optimal_h = " || glob_reached_optimal_h || " iter = " || iter || " glob_h = " || glob_h ); if glob_reached_optimal_h = 0 then est_err := estimate_trunc_err(); print("Estimated truncation error = " , est_err); print("Maximum truncation error = " , max_trunc_err); if (est_err < max_trunc_err) and (hnew <> glob_hmax) then if iter = 0 then glob_unchanged_h_cnt := 0; glob_reached_optimal_h := 0; hnew := 2.0 * glob_h; if hnew > glob_hmax then hnew := glob_hmax; fi; glob_last_good_h := hnew; glob_correct_start_flag := 1; print("Solution starting with Initally determined glob_h = " , hnew); cs_info("Initial increase h = " , hnew ); break; else if iter > 0 then print("Restarting to try larger glob_h!"); cs_info("Restart h = " , hnew); hnew := 2.0 * glob_h; if hnew > glob_hmax then hnew := glob_hmax; fi; glob_reached_optimal_h := 0; glob_last_good_h := hnew; glob_correct_start_flag := 1; break; fi; #revised 8/8/8 fi; # Guess 1 fi; # Guess 2 fi; # Guess 3 if (est_err < max_trunc_err) then print("glob_h unchanged"); glob_last_good_h := glob_h; glob_unchanged_h_cnt = glob_unchanged_h_cnt + 1; if (glob_unchanged_h_cnt >= MAX_UNCHANGED) then print("Reached Optimal"); glob_reached_optimal_h := 1; glob_curr_iter_when_opt := current_iter; glob_correct_start_flag := 0; glob_optimal_clock_start_sec := time(); glob_optimal_start := glob_x[1]; else print("Backtracking to previous glob_h!"); print("Reached Optimal"); hnew = glob_h / 2.0; cs_info("Backtrack h = " , hnew); glob_reached_optimal_h := 1; glob_curr_iter_when_opt := current_iter; glob_optimal_clock_start_sec := time(); glob_optimal_start := glob_x_start; glob_correct_start_flag := 1; break; fi; fi; fi; glob_x[1] := glob_x[1]+(glob_h); glob_x[2] := (glob_h); jump_glob_y(); if ((iter mod glob_progress_ind) = 0) then print_alot(1,glob_x,glob_y_higher,3,"glob_y","glob_x"); fi; glob_hsav := (glob_h); hnewsav := (hnew); if hnew < glob_h and hnew > glob_hmin_init then hnew := (adjust_step(glob_x,glob_y_higher,glob_hsav,hnewsav,3)); fi; cs_info("Before potential break iter == " , iter ); glob_clock_sec := time(); glob_h := hnew; cs_info("before end cs loop adj"); if glob_reached_optimal_h = 1 or iter >= glob_max_iter or (time() - glob_orig_start_sec) >= glob_max_sec then glob_correct_start_flag = 0; fi; cs_info("after end cs loop adj iter = " , iter); # djd_cnt := djd_cnt + 1; # if djd_cnt > 30 then break; fi; od; # Guess 3 # djd_cnt := djd_cnt + 1; # if djd_cnt > 30 then break; fi; od; # Guess 4 cs_info("After All loops"); print("Finished!"); fprintf(glob_out_fd,"Finished!\n"); if iter >= glob_max_iter then print("Maximum Iterations Reached!"); fprintf(glob_out_fd,"Maximum Iterations Reached!\n"); fi; if (time() - glob_orig_start_sec) >= glob_max_sec then print("Maximum Time Reached!"); fprintf(glob_out_fd,"Maximum Time Reached!\n"); fi; print("Output in diff2.ode.out.txt"); fclose(glob_progress_fd); fclose(glob_out_fd); # END OUTFILEMAIN