[Tutor] hi
Vick
vick1975 at orange.mu
Sun Aug 11 07:09:43 CEST 2013
> -----Original Message-----
> From: Oscar Benjamin [mailto:oscar.j.benjamin at gmail.com]
> Sent: Saturday, 10 August, 2013 19:10
> To: Vick
> Cc: tutor at python.org
> Subject: Re: [Tutor] hi
>
>
> How long a time it takes depends on many different things: the language,
the
> precision, how good a programmer you are. Much of the "this integrator is
> better than that one" arguments are implicitly tied to certain
presumptions
> about how your code gets implemented.
[Vick] Well, if you coded the methods in the same format, structure and
language, surely the higher order method will produce better accuracy for
the same ODE problem and for the same amount of time than the lower order.
> Here is an example that shows we can achieve the precision of scipy's
> integrator (for the system you posted) easily using a standard 4th order
> Runge-Kutta:
[Vick] Your code shows just that. I ran it and it went on for like about 40
minutes. In that amount of time, the memory was taxed a lot. For the same
ODE problem, RK4 should produce an error of 1e-4 in a reasonable amount of
time (about 2 secs) and utilizing normal computer memory for its
calculations. As I mentioned before I used this problem as a test for
accuracy to compare different methods. I use percentage change (x 100) to
derive the error for each method. RK4 in this instance scores about 1e-4
accuracy. With this problem I've tested RKF5, DPRK45, some 6 and 7 order RK
methods and the last one is the DOPRI8(7)13.
The Leapfrog method is a symplectic integrator and I am currently trying to
code it. I have it in Excel format already.
The following is my code for the ODE problem test for RK4 and DOPRI8(7)13:
(Note that they are using the same language, structure and format) (you will
also need mpmath to run: I use mpmath for large number calculations)
The code should run using less than 2 secs and the output is pre-defined for
time, y , true value and error. The error for the two methods will show the
difference between the different orders. You should change the call to each
method for comparison as I coded it to run one method at any one time.
from mpmath import *
mp.dps=30
def drange(start, stop, step):
r = start
while r < stop:
yield r
r += step
def f(t, y):
return (mpf('4') * (exp(1) ** (mpf('0.8') * t))) - (mpf('0.5') * y)
def rungk4(t, y, dt):
k1 = dt * f(t, y)
k2 = dt * f(t + dt / mpf('2'), y + k1 / mpf('2'))
k3 = dt * f(t + dt / mpf('2'), y + k2 / mpf('2'))
k4 = dt * f(t + dt, y + k3)
return y + (k1 + mpf('2') * (k2 + k3) + k4) / mpf('6')
def rungkdp8(t, y, dt):
k1 = dt * f(t, y)
k2 = dt * f(t + dt / mpf('18'), y + k1 / mpf('18'))
k3 = dt * f(t + dt / mpf('12'), y + k1 / mpf('48') + k2 / mpf('16'))
k4 = dt * f(t + dt / mpf('8'), y + k1 / mpf('32') + k3 * mpf('3') /
mpf('32'))
k5 = dt * f(t + dt * mpf('5') / mpf('16'), y + k1 * mpf('5') / mpf('16')
- k3 * mpf('75') / mpf('64') + k4 * mpf('75') / mpf('64'))
k6 = dt * f(t + dt * mpf('3') / mpf('8'), y + k1 * mpf('3') / mpf('80')
+ k4 * mpf('3') / mpf('16') + k5 * mpf('3') / mpf('20'))
k7 = dt * f(t + dt * mpf('59') / mpf('400'), y + k1 * mpf('29443841') /
mpf('614563906') + k4 * mpf('77736538') / mpf('692538347') - k5 *
mpf('28693883') / mpf('1125000000') + k6 * mpf('23124283') /
mpf('1800000000'))
k8 = dt * f(t + dt * mpf('93') / mpf('200'), y + k1 * mpf('16016141') /
mpf('946692911') + k4 * mpf('61564180') / mpf('158732637') + k5 *
mpf('22789713') / mpf('633445777') + k6 * mpf('545815736') /
mpf('2771057229') - k7 * mpf('180193667') / mpf('1043307555'))
k9 = dt * f(t + dt * mpf('5490023248') / mpf('9719169821'), y + k1 *
mpf('39632708') / mpf('573591083') - k4 * mpf('433636366') /
mpf('683701615') - k5 * mpf('421739975') / mpf('2616292301') + k6 *
mpf('100302831') / mpf('723423059') + k7 * mpf('790204164') /
mpf('839813087') + k8 * mpf('800635310') / mpf('3783071287'))
k10 = dt * f(t + dt * mpf('13') / mpf('20'), y + k1 * mpf('246121993') /
mpf('1340847787') - k4 * mpf('37695042795') / mpf('15268766246') - k5 *
mpf('309121744') / mpf('1061227803') - k6 * mpf('12992083') /
mpf('490766935') + k7 * mpf('6005943493') / mpf('2108947869') + k8 *
mpf('393006217') / mpf('1396673457') + k9 * mpf('123872331') /
mpf('1001029789'))
k11 = dt * f(t + dt * mpf('1201146811') / mpf('1299019798'), y - k1 *
mpf('1028468189') / mpf('846180014') + k4 * mpf('8478235783') /
mpf('508512852') + k5 * mpf('1311729495') / mpf('1432422823') - k6 *
mpf('10304129995') / mpf('1701304382') - k7 * mpf('48777925059') /
mpf('3047939560') + k8 * mpf('15336726248') / mpf('1032824649') - k9 *
mpf('45442868181') / mpf('3398467696') + k10 * mpf('3065993473') /
mpf('597172653'))
k12 = dt * f(t + dt, y + k1 * mpf('185892177') / mpf('718116043') - k4 *
mpf('3185094517') / mpf('667107341') - k5 * mpf('477755414') /
mpf('1098053517') - k6 * mpf('703635378') / mpf('230739211') + k7 *
mpf('5731566787') / mpf('1027545527') + k8 * mpf('5232866602') /
mpf('850066563') - k9 * mpf('4093664535') / mpf('808688257') + k10 *
mpf('3962137247') / mpf('1805957418') + k11 * mpf('65686358') /
mpf('487910083'))
k13 = dt * f(t + dt, y + k1 * mpf('403863854') / mpf('491063109') - k4 *
mpf('5068492393') / mpf('434740067') - k5 * mpf('411421997') /
mpf('543043805') + k6 * mpf('652783627') / mpf('914296604') + k7 *
mpf('11173962825') / mpf('925320556') - k8 * mpf('13158990841') /
mpf('6184727034') + k9 * mpf('3936647629') / mpf('1978049680') - k10 *
mpf('160528059') / mpf('685178525') + k11 * mpf('248638103') /
mpf('1413531060'))
return y + (k1 * mpf('14005451') / mpf('335480064') - k6 *
mpf('59238493') / mpf('1068277825') + k7 * mpf('181606767') /
mpf('758867731') + k8 * mpf('561292985') / mpf('797845732') - k9 *
mpf('1041891430') / mpf('1371343529') + k10 * mpf('760417239') /
mpf('1151165299') + k11 * mpf('118820643') / mpf('751138087') - k12 *
mpf('528747749') / mpf('2220607170') + k13 / mpf('4'))
#return y + (k1 * mpf('13451932') / mpf('455176623') - k6 *
mpf('808719846') / mpf('976000145') + k7 * mpf('1757004468') /
mpf('5645159321') + k8 * mpf('656045339') / mpf('265891186') - k9 *
mpf('3867574721') / mpf('1518517206') + k10 * mpf('465885868') /
mpf('322736535') + k11 * mpf('53011238') / mpf('667516719') + k12 * mpf('2')
/ mpf('45'))
h = mpf('0.25')
t = mpf('0')
y = mpf('2')
tn = mpf('4.25')
tv
=(mpf('4')/mpf('1.3'))*exp(mpf('0.8')*t)-(mpf('1.4')/mpf('1.3'))*exp(mpf('-0
.5')*t)
ti = t+h
tv1
=(mpf('4')/mpf('1.3'))*exp(mpf('0.8')*ti)-(mpf('1.4')/mpf('1.3'))*exp(mpf('-
0.5')*ti)
mi = rungkdp8(t,y,h) #change method here (either rungkdp8 or rungk4)
er1 =abs((tv1-mi)/tv1)*mpf('100')
print "x Y True Value
Error"
print
"---------------------------------------------------------------------------
----------------------"
print t," ",y," ",tv
print ti," ", mi," ",tv1," ",er1
for i in drange (ti, tn,h):
mi = rungkdp8(i,mi,h) #change method here (either rungkdp8 or rungk4)
tvn
=(mpf('4')/mpf('1.3'))*exp(mpf('0.8')*(i+h))-(mpf('1.4')/mpf('1.3'))*exp(mpf
('-0.5')*(i+h))
er2 = abs((tvn-mi)/tvn)*mpf('100')
print i+h," ", mi," ",tvn," ",er2
Vick
More information about the Tutor
mailing list