Changeset 57eb6a4 in sasmodels


Ignore:
Timestamp:
Oct 17, 2017 4:29:59 PM (4 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
2a602c7
Parents:
0d19f42
Message:

add hints for generating taylor series code from sympy

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/precision.py

    r487e695 r57eb6a4  
    11#!/usr/bin/env python 
    22r""" 
    3 Show numerical precision of $2 J_1(x)/x$. 
     3Show numerical precision of various expressions. 
     4 
     5Evaluates the same function(s) in single and double precision and compares 
     6the results to 500 digit mpmath evaluation of the same function. 
     7 
     8Note: a quick way to generation C and python code for taylor series 
     9expansions from sympy: 
     10 
     11    import sympy as sp 
     12    x = sp.var("x") 
     13    f = sp.sin(x)/x 
     14    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term 
     15    p = sp.horner(t)   # Horner representation 
     16    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero 
     17    p = p.n(15)  # evaluate coefficients to 15 digits (optional) 
     18    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code 
     19    py_code = c[:-1]  # strip semicolon to convert c to python 
     20 
     21    # mpmath has pade() rational function approximation, which might work 
     22    # better than the taylor series for some functions: 
     23    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M) 
     24    P = sum(a*x**n for n,a in enumerate(reversed(P))) 
     25    Q = sum(a*x**n for n,a in enumerate(reversed(Q))) 
     26    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p")) 
     27 
     28    # There are richardson and shanks series accelerators in both sympy 
     29    # and mpmath that may be helpful. 
    430""" 
    531from __future__ import division, print_function 
Note: See TracChangeset for help on using the changeset viewer.