Hi Doug

 

thank for your post.

But how'd you do, for exemple, if you were in one or more of these situations (witouth using for cycles)?

1)  You want to check, before the assignment, wether 3^(i-1) (in your ex) is positive or negative, or smaller than 1, etc, and print a string in a file that takes memory of this checks

2) Suppose that your assigment function is composed of many addends, and you want, before the assigment is done for a particular set of indexes, to check for the positiveness of all the addends, and then add all the positives togher, all the negatives together, and finally this two quantities togher

3) You want to assign only one part of the entries, eg only these for whom the second index is smaller that the first.

Maybe I'm wrong, but I've always thought that the direct assigment is possible only in very simple cases...
Salvo

Ok, the problem should really be the Integration. let me try to make a little resume: the main point of my program is a scalar product IntegrationTools:-Expand( 2.0*Int(f^a/S_h(f), f=20.. f_max ))); with S_h given before and a on the range [-8..8] with 1/3 steps (-8, -23/3 ... 23/3, +8) f_max is not constant, but depends on the mass of the system I'm considering, in this way: f_max= evalf(6^(3/2)*Pi*M_tot)^(-1): where M_tot is in units of solar mass M_solar:=4.927*10^(-6), so M_tot = (a positive number * M_solar). My program begins to give strange or unstable results for big M_tot (ie: small f_high). here for ex the same quantity calculated 200 times for different values of M_tot, from 10*M_solar to 130*M_solar 10: 30: 60: 90: 130: so now I'm pretty sure that the problem is indeed the integration with 15 digits, that becomes more and more inaccurate when the up limits of integration approaches the lower one. I'd be nice to understand why there are these "plateau" of some dozens of point, where the results are quite the same, before the next jump. I'd stress that these are only intermediary results, that I'll manipulate later in my code, this explain, I think, why mi final results are so unstable. Suggestions welcomed!!! Tomorrow I'll try the methods Alex proposed before. Thanks all Salvo

Ok, the problem should really be the Integration. let me try to make a little resume: the main point of my program is a scalar product IntegrationTools:-Expand( 2.0*Int(f^a/S_h(f), f=20.. f_max ))); with S_h given before and a on the range [-8..8] with 1/3 steps (-8, -23/3 ... 23/3, +8) f_max is not constant, but depends on the mass of the system I'm considering, in this way: f_max= evalf(6^(3/2)*Pi*M_tot)^(-1): where M_tot is in units of solar mass M_solar:=4.927*10^(-6), so M_tot = (a positive number * M_solar). My program begins to give strange or unstable results for big M_tot (ie: small f_high). here for ex the same quantity calculated 200 times for different values of M_tot, from 10*M_solar to 130*M_solar 10: 30: 60: 90: 130: so now I'm pretty sure that the problem is indeed the integration with 15 digits, that becomes more and more inaccurate when the up limits of integration approaches the lower one. I'd be nice to understand why there are these "plateau" of some dozens of point, where the results are quite the same, before the next jump. I'd stress that these are only intermediary results, that I'll manipulate later in my code, this explain, I think, why mi final results are so unstable. Suggestions welcomed!!! Tomorrow I'll try the methods Alex proposed before. Thanks all Salvo

Here there is a plot of the 200 values of the variable before (.94 etc), I can't see any visible trend that would stand the bug code hypothesis...

Here there is a plot of the 200 values of the variable before (.94 etc), I can't see any visible trend that would stand the bug code hypothesis...

Hi Alex, thank for you rapidity. Actually, I don't need 30 digits in my final result (that would be crazy, and usefulness for my purposes), I'll be happy of having result with 3-4 "true and faithful" digits. The problem is that I get stuff like this sequence (the same quantity in different runs of my code): .9401997463890721916951909563724e-1/rho .9472918729980896437917381518568e-1/rho .9500773424747221984237661836773e-1/rho .9472950697941090283626108434313e-1/rho .9472940828757043301072377011556e-1/rho .9472961281854383076949160772922e-1/rho .9392101409109034551742078411885e-1/rho .9472943656328361000186837799826e-1/rho .9521506980238236412044032417133e-1/rho .9472951048695593323602732377927e-1/rho .9636560726146277227732486502852e-1/rho .9472944126212900040706363608476e-1/rho .9472958244904912577474031472477e-1/rho etc. You can see that sometime the second digits is unstable. I'd like to understand if the problem come from the integral, which is the only part of my code which still works with 15 digits because of the time, or if I've made a big mistake somewhere. I don't know if I can copy here a part of my code, that I've changed following your suggestion (proc declaration before, calling inside the for), it's long and rather complicated to follow...

Hi Alex, thank for you rapidity. Actually, I don't need 30 digits in my final result (that would be crazy, and usefulness for my purposes), I'll be happy of having result with 3-4 "true and faithful" digits. The problem is that I get stuff like this sequence (the same quantity in different runs of my code): .9401997463890721916951909563724e-1/rho .9472918729980896437917381518568e-1/rho .9500773424747221984237661836773e-1/rho .9472950697941090283626108434313e-1/rho .9472940828757043301072377011556e-1/rho .9472961281854383076949160772922e-1/rho .9392101409109034551742078411885e-1/rho .9472943656328361000186837799826e-1/rho .9521506980238236412044032417133e-1/rho .9472951048695593323602732377927e-1/rho .9636560726146277227732486502852e-1/rho .9472944126212900040706363608476e-1/rho .9472958244904912577474031472477e-1/rho etc. You can see that sometime the second digits is unstable. I'd like to understand if the problem come from the integral, which is the only part of my code which still works with 15 digits because of the time, or if I've made a big mistake somewhere. I don't know if I can copy here a part of my code, that I've changed following your suggestion (proc declaration before, calling inside the for), it's long and rather complicated to follow...

Hi Alex there is something I don't understand in your post. You say that I should use Digits:=evalhf(Digits) ie Digits:=15, in the machines I run the code. But this is exactly what I've done before opening this thread, and I had good results (for not too small upper integration extrema). What I want is exactly to use Digits>=16 without loosing all these performances. Today I've tried this My_Int:= proc(func, f_min,f_max) Digits:=16; bins:=100; tay_order:=10; delta_f:= (f_max-f_min)/bins; f_low:=f_min; f_up:=f_min+delta_f; temp:=0; tot:=0; while f_up <=f_max do middle_f:= f_low+delta_f/2; tayl:= convert(series(S_h(f), f=middle_f, tay_order),polynom); temp:=evalf(Int( func/S_h(f), f=f_low..f_up); tot:=tot+temp; f_low:=f_up; f_up:=f_up+delta_f; end do; return tot; end proc; where eventually the values of bins and tay_order can be raised depending on the amplitude of f_max-f_min. It gives good results, and is relatively fast (much faster that evalf[16](Int(func/S_h, f=f_min..f_max)) Otherwise, is there in the Matlab package of maple the numerical integration of matlab (quad, or quadl)? I've passed my functions to a friend and he verified that matlab if much faster for this particular integral, with 32 digits. Thanks S.

Hi Alex there is something I don't understand in your post. You say that I should use Digits:=evalhf(Digits) ie Digits:=15, in the machines I run the code. But this is exactly what I've done before opening this thread, and I had good results (for not too small upper integration extrema). What I want is exactly to use Digits>=16 without loosing all these performances. Today I've tried this My_Int:= proc(func, f_min,f_max) Digits:=16; bins:=100; tay_order:=10; delta_f:= (f_max-f_min)/bins; f_low:=f_min; f_up:=f_min+delta_f; temp:=0; tot:=0; while f_up <=f_max do middle_f:= f_low+delta_f/2; tayl:= convert(series(S_h(f), f=middle_f, tay_order),polynom); temp:=evalf(Int( func/S_h(f), f=f_low..f_up); tot:=tot+temp; f_low:=f_up; f_up:=f_up+delta_f; end do; return tot; end proc; where eventually the values of bins and tay_order can be raised depending on the amplitude of f_max-f_min. It gives good results, and is relatively fast (much faster that evalf[16](Int(func/S_h, f=f_min..f_max)) Otherwise, is there in the Matlab package of maple the numerical integration of matlab (quad, or quadl)? I've passed my functions to a friend and he verified that matlab if much faster for this particular integral, with 32 digits. Thanks S.

Hi Axel if I plot f^a/S_h(f) with a on a range [-3 .. 3] (this is the range of my interest) and with f in [20 .. 2000] I see that there are not infinities. The only zero of the denominator should be for 0