math_.sympy_with_test のソースコード

import argparse
import numpy as np
import sympy as sym
from sympy.parsing.sympy_parser import parse_expr
import sys


"""
python sympy_with_test.py "1/(exp(x)+1)" --mode diff --norder 8 --xmin -2 --xmax 2 --nx 5
python sympy_with_test.py "exp(x)" --mode taylor --norder 8 --xmin -2 --xmax 2 --nx 5
python sympy_with_test.py "x**2 + cos(x)" --mode int --export
"""

import argparse
import numpy as np
import sympy as sym
from sympy.parsing.sympy_parser import parse_expr
import sys



def get_calculation_functions(formula_str, mode, x0, norder):
    # ライブラリ名なしの関数対応テーブル
    func_map = {
        'sin': np.sin, 'cos': np.cos, 'tan': np.tan,
        'exp': np.exp, 'log': np.log, 'sqrt': np.sqrt,
        'abs': np.abs, 'sinh': np.sinh, 'cosh': np.cosh
    }

    x = sym.symbols('x')
    try:
        expr = parse_expr(formula_str)
        
        if mode == "diff":
            target_expr = sym.diff(expr, x)
            label = "Differentiation (1st)"
        elif mode == "int":
            target_expr = sym.integrate(expr, x)
            label = "Integration (Indefinite)"
        elif mode == "taylor":
            target_expr = sym.series(expr, x, x0, norder).removeO()
            label = f"Taylor Expansion (n={norder}, x0={x0})"
        
        # 数値計算用関数の生成
        f_num = sym.lambdify(x, expr, modules=[func_map, "numpy"])
        target_num = sym.lambdify(x, target_expr, modules=[func_map, "numpy"])
        
        return expr, target_expr, f_num, target_num, label
    except Exception as e:
        print(f"解析エラー: {e}")
        sys.exit(1)




def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("formula", type=str)
    parser.add_argument("--mode", choices=["diff", "int", "taylor"], default="diff")
    parser.add_argument("--xmin", type=float, default=-1.0)
    parser.add_argument("--xmax", type=float, default=1.0)
    parser.add_argument("--nx", type=int, default=11)
    parser.add_argument("--norder", type=int, default=5)
    parser.add_argument("--x0", type=float, default=0.0)
    parser.add_argument("-s", "--step", type=float, default=1e-5)
    parser.add_argument("--export", action="store_true")

    args = parser.parse_args()
    expr, t_expr, f_func, t_func, label = get_calculation_functions(
        args.formula, args.mode, args.x0, args.norder
    )

    h = args.step
    x_vals = np.linspace(args.xmin, args.xmax, args.nx)

    print(f"\n{'='*100}")
    print(f" Mode: {label}")
    print(f" Result Formula: {t_expr}")
    print(f"{'='*100}")
    
    # モードに応じたヘッダー
    if args.mode == "diff":
        col_name = "Numeric Diff"
    elif args.mode == "int":
        col_name = "Numeric Integ*" # 簡易的な確認
    else:
        col_name = "Original f(x)"

    print(f"{'x':>8} | {'Original f(x)':>18} | {col_name:>18} | {'SymPy Result':>18} | {'Rel.Err':>10}")
    print("-" * 100)

    for xv in x_vals:
        try:
            v_orig = f_func(xv)
            v_sympy = t_func(xv)
            
            # 比較用の数値計算(Transformed Value)
            if args.mode == "diff":
                # 中央差分
                v_trans = (f_func(xv + h) - f_func(xv - h)) / (2 * h)
            elif args.mode == "int":
                # 不定積分の比較は難しいが、ここでは x0 からの定積分(台形則)で近似値を出す例
                # (簡易的に微分して元の関数に戻るか見る方が正確ですが、要望に合わせ数値比較)
                v_trans = float('nan') # 積分は評価が特殊なため一旦nan
            elif args.mode == "taylor":
                v_trans = v_orig # テイラー展開は元の関数値と比較

            # 誤差計算 (SymPyの結果 vs 数値的アプローチ)
            err_ref = v_trans if not np.isnan(v_trans) else v_orig
            rel_err = abs(v_sympy - err_ref) / abs(err_ref) if err_ref != 0 else 0

            print(f"{xv:8.3f} | {v_orig:18.8f} | {v_trans:18.8f} | {v_sympy:18.8f} | {rel_err:10.2e}")
        except:
            print(f"{xv:8.3f} | 計算エラー")

    if args.export:
        # 前述のコード出力関数をここに呼ぶ
        pass


if __name__ == "__main__":
    main()