Re: question for lyuka and the Ostrowski method  bumped! Message #4 Posted by Lyuka on 25 Aug 2011, 10:01 p.m., in response to message #1 by hugh steers
Hi hugh steers, Namir and all who interested.
I've modified the code (Rev.2.94) to struggle to escape from the flat region or extremes, this morning.
Now, it can solve the water saturation pressure equation started at the very flat region e.g. [273.15, 273.0].
Please try the latest code (Rev.2.94).
Quote:
/*s(x) : Water saturation pressures (Sonntag)*/
double
s(x)
double x;
{
double t;
double p = 1e5;
t = x + 273.15;
if (t < 0) return p;
return exp(((16.635794 + 2.433502 * log(t) + 0.00001673952 * t * t)  0.02711193 * t)
 6096.9385 / t)  p;
}
main: s(t) is 'P_water_sat_Sonntag'
ost:a=273.15, b=273.1, c=273.125, f=1.000000000000e05, i=0, p=0
ost:a=273.15, b=273.225, c=273.206078387, f=1.000000000000e05, i=1, p=1
ost:a=273.15, b=273.0375, c=273.064372309, f=1.000000000000e05, i=2, p=2
ost:a=273.15, b=273.31875, c=273.295157305, f=1.000000000000e05, i=3, p=3
ost:a=273.15, b=272.896875, c=273.090425146, f=1.000000000000e05, i=4, p=4
ost:a=273.15, b=273.5296875, c=273.285100373, f=1.000000000000e05, i=5, p=5
ost:a=273.15, b=272.58046875, c=272.731564881, f=1.000000000000e05, i=6, p=6
ost:a=273.15, b=274.004296875, c=273.411032864, f=1.000000000000e05, i=7, p=7
ost:a=273.15, b=271.868554688, c=272.448762741, f=1.000000000000e05, i=8, p=8
ost:a=273.15, b=275.072167969, c=274.073068154, f=1.000000000000e05, i=9, p=9
ost:a=273.15, b=270.266748047, c=271.383252597, f=1.000000000000e05, i=10, p=10
ost:a=273.15, b=277.47487793, c=274.800747146, f=1.000000000000e05, i=11, p=11
ost:a=273.15, b=266.662683105, c=269.830272597, f=1.000000000000e05, i=12, p=12
ost:a=273.15, b=282.880975342, c=281.866044434, f=1.000000000000e05, i=13, p=13
ost:a=273.15, b=258.553536987, c=260.536164572, f=1.000000000000e05, i=14, p=14
ost:a=273.15, b=295.044694519, c=286.697293227, f=1.000000000000e05, i=15, p=15
ost:a=273.15, b=240.307958221, c=250.48456364, f=1.000000000000e05, i=16, p=16
ost:a=273.15, b=322.413062668, c=282.332745028, f=1.000000000000e05, i=17, p=17
ost:a=273.15, b=199.255405998, c=229.286334046, f=1.000000000000e05, i=18, p=18
ost:a=273.15, b=383.991891003, c=281.658556097, f=1.000000000000e05, i=19, p=19
ost:a=273.15, b=106.887163496, c=227.42340973, f=1.000000000000e05, i=20, p=20
ost:a=273.15, b=106.887163496, c=242.794125415, f=1.000000000000e05, i=21, p=21
ost:a=273.15, b=106.887163496, c=145.766953091, f=9.999848970363e06, i=22, p=22
ost:a=106.887163496, b=145.766953091, c=106.88706331, f=1.155857857677e06, i=23, p=22
ost:a=145.766953091, b=106.88706331, c=106.338424292, f=6.002266450796e08, i=24, p=22
ost:a=106.88706331, b=106.338424292, c=106.311453476, f=2.988926674806e09, i=25, p=22
ost:a=106.338424292, b=106.311453476, c=106.310044116, f=1.842629348946e13, i=26, p=22
ost:a=106.311453476, b=106.310044116, c=106.310044029, f=1.496866036476e20, i=27, p=22
main: 273.15, 273.1, s(106.310044029) = 1.49686603648e20, r=1
And also, it can find a root of (x1)^5 in 27 iterations.
Quote:
/* f6(x) = (x1)^5 */
ost:a=0, b=1.5, c=0.75, f=9.765625000000e04, i=0, p=0
ost:a=1.5, b=0.75, c=0.761904761905, f=7.651622719419e04, i=1, p=0
ost:a=0.75, b=0.761904761905, c=0.803189200351, f=2.952872052788e04, i=2, p=0
ost:a=0.761904761905, b=0.803189200351, c=0.850693877302, f=7.419729755354e05, i=3, p=0
ost:a=0.803189200351, b=0.850693877302, c=0.878293985111, f=2.670300955364e05, i=4, p=0
ost:a=0.850693877302, b=0.878293985111, c=0.902610176346, f=8.761286801440e06, i=5, p=0
ost:a=0.878293985111, b=0.902610176346, c=0.92234070912, f=2.824666003999e06, i=6, p=0
ost:a=0.902610176346, b=0.92234070912, c=0.937768603186, f=9.333569060579e07, i=7, p=0
ost:a=0.92234070912, b=0.937768603186, c=0.950216234015, f=3.058008812555e07, i=8, p=0
ost:a=0.937768603186, b=0.950216234015, c=0.960176299828, f=1.001631669682e07, i=9, p=0
ost:a=0.950216234015, b=0.960176299828, c=0.968133427318, f=3.286069794729e08, i=10, p=0
ost:a=0.960176299828, b=0.968133427318, c=0.974504259003, f=1.077303803554e08, i=11, p=0
ost:a=0.968133427318, b=0.974504259003, c=0.979601206998, f=3.532013499117e09, i=12, p=0
ost:a=0.974504259003, b=0.979601206998, c=0.983678883914, f=1.158108552378e09, i=13, p=0
ost:a=0.979601206998, b=0.983678883914, c=0.986941586119, f=3.797100974770e10, i=14, p=0
ost:a=0.983678883914, b=0.986941586119, c=0.989552033388, f=1.244969983086e10, i=15, p=0
ost:a=0.986941586119, b=0.989552033388, c=0.991640629072, f=4.081952835848e11, i=16, p=0
ost:a=0.989552033388, b=0.991640629072, c=0.993311710543, f=1.338367252555e11, i=17, p=0
ost:a=0.991640629072, b=0.993311710543, c=0.994648733138, f=4.388165984146e12, i=18, p=0
ost:a=0.993311710543, b=0.994648733138, c=0.995718477996, f=1.438768435661e12, i=19, p=0
ost:a=0.994648733138, b=0.995718477996, c=0.996574375736, f=4.717355769657e13, i=20, p=0
ost:a=0.995718477996, b=0.996574375736, c=0.99725917517, f=1.546701119566e13, i=21, p=0
ost:a=0.996574375736, b=0.99725917517, c=0.997807079768, f=5.071240091293e14, i=22, p=0
ost:a=0.99725917517, b=0.997807079768, c=0.998245455502, f=1.662730778483e14, i=23, p=0
ost:a=0.997807079768, b=0.998245455502, c=0.998596197729, f=5.451671807322e15, i=24, p=0
ost:a=0.998245455502, b=0.998596197729, c=0.99887682483, f=1.787464685796e15, i=25, p=0
ost:a=0.998596197729, b=0.99887682483, c=0.999101353168, f=5.860642591876e16, i=26, p=0
ost:a=0.99887682483, b=0.999101353168, c=0.999280997168, f=1.921555815013e16, i=27, p=0
main:f6: 0, 1.5, s(0.999280997168) = 1.92155581501e16, r=1
Regards,
Lyuka
Edited: 25 Aug 2011, 11:59 p.m.
