<s> from bisect import bisect <EOL> from . . libmp . backend import xrange <EOL> class ODEMethods ( object ) : <EOL> pass <EOL> def ode_taylor ( ctx , derivs , x0 , y0 , tol_prec , n ) : <EOL> h = tol = ctx . ldexp ( <NUM_LIT:1> , - tol_prec ) <EOL> dim = len ( y0 ) <EOL> xs = [ x0 ] <EOL> ys = [ y0 ] <EOL> x = x0 <EOL> y = y0 <EOL> orig = ctx . prec <EOL> try : <EOL> ctx . prec = orig * ( <NUM_LIT:1> + n ) <EOL> for i in range ( n ) : <EOL> fxy = derivs ( x , y ) <EOL> y = [ y [ i ] + h * fxy [ i ] for i in xrange ( len ( y ) ) ] <EOL> x += h <EOL> xs . append ( x ) <EOL> ys . append ( y ) <EOL> ser = [ [ ] for d in range ( dim ) ] <EOL> for j in range ( n + <NUM_LIT:1> ) : <EOL> s = [ <NUM_LIT:0> ] * dim <EOL> b = ( - <NUM_LIT:1> ) ** ( j & <NUM_LIT:1> ) <EOL> k = <NUM_LIT:1> <EOL> for i in range ( j + <NUM_LIT:1> ) : <EOL> for d in range ( dim ) : <EOL> s [ d ] += b * ys [ i ] [ d ] <EOL> b = ( b * ( j - k + <NUM_LIT:1> ) ) // ( - k ) <EOL> k += <NUM_LIT:1> <EOL> scale = h ** ( - j ) / ctx . fac ( j ) <EOL> for d in range ( dim ) : <EOL> s [ d ] = s [ d ] * scale <EOL> ser [ d ] . append ( s [ d ] ) <EOL> finally : <EOL> ctx . prec = orig <EOL> radius = ctx . one <EOL> for ts in ser : <EOL> if ts [ - <NUM_LIT:1> ] : <EOL> radius = min ( radius , ctx . nthroot ( tol / abs ( ts [ - <NUM_LIT:1> ] ) , n ) ) <EOL> radius /= <NUM_LIT:2> <EOL> return ser , x0 + radius <EOL> def odefun ( ctx , F , x0 , y0 , tol = None , degree = None , method = '<STR_LIT>' , verbose = False ) : <EOL> r"""<STR_LIT>""" <EOL> if tol : <EOL> tol_prec = int ( - ctx . log ( tol , <NUM_LIT:2> ) ) + <NUM_LIT:10> <EOL> else : <EOL> tol_prec = ctx . prec + <NUM_LIT:10> <EOL> degree = degree or ( <NUM_LIT:3> + int ( <NUM_LIT:3> * ctx . dps / <NUM_LIT> ) ) <EOL> workprec = ctx . prec + <NUM_LIT> <EOL> try : <EOL> len ( y0 ) <EOL> return_vector = True <EOL> except TypeError : <EOL> F_ = F <EOL> F = lambda x , y : [ F_ ( x , y [ <NUM_LIT:0> ] ) ] <EOL> y0 = [ y0 ] <EOL> return_vector = False <EOL> ser , xb = ode_taylor ( ctx , F , x0 , y0 , tol_prec , degree ) <EOL> series_boundaries = [ x0 , xb ] <EOL> series_data = [ ( ser , x0 , xb ) ] <EOL> def mpolyval ( ser , a ) : <EOL> return [ ctx . polyval ( s [ : : - <NUM_LIT:1> ] , a ) for s in ser ] <EOL> def get_series ( x ) : <EOL> if x < x0 : <EOL> raise ValueError <EOL> n = bisect ( series_boundaries , x ) <EOL> if n < len ( series_boundaries ) : <EOL> return series_data [ n - <NUM_LIT:1> ] <EOL> while <NUM_LIT:1> : <EOL> ser , xa , xb = series_data [ - <NUM_LIT:1> ] <EOL> if verbose : <EOL> print ( "<STR_LIT>" % ( xa , xb ) ) <EOL> y = mpolyval ( ser , xb - xa ) <EOL> xa = xb <EOL> ser , xb = ode_taylor ( ctx , F , xb , y , tol_prec , degree ) <EOL> series_boundaries . append ( xb ) <EOL> series_data . append ( ( ser , xa , xb ) ) <EOL> if x <= xb : <EOL> return series_data [ - <NUM_LIT:1> ] <EOL> def interpolant ( x ) : <EOL> x = ctx . convert ( x ) <EOL> orig = ctx . prec <EOL> try : <EOL> ctx . prec = workprec <EOL> ser , xa , xb = get_series ( x ) <EOL> y = mpolyval ( ser , x - xa ) <EOL> finally : <EOL> ctx . prec = orig <EOL> if return_vector : <EOL> return [ + yk for yk in y ] <EOL> else : <EOL> return + y [ <NUM_LIT:0> ] <EOL> return interpolant <EOL> ODEMethods . odefun = odefun <EOL> if __name__ == "<STR_LIT:__main__>" : <EOL> import doctest <EOL> doctest . testmod ( ) </s>
