pythonsympyleast-squaresmpmathfunction-approximation

Using Python to find best square approximation of Piecewise Functions？

Given a parabola, $f(x)=10(x-1)^2-1,x\in \Omega =[1,2]$, find the best approximation in the space of all linear functions. The Python code is as follows:

import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
from mpmath import *

def least_squares(f, psi, Omega, symbolic=True):
    N = len(psi) - 1
    A = sym.zeros(N+1, N+1)
    b = sym.zeros(N+1, 1)
    x = sym.Symbol('x')
    for i in range(N+1):
        for j in range(i, N+1):
            integrand = psi[i]*psi[j]
            if symbolic:
                I = sym.integrate(integrand, (x, Omega[0], Omega[1]))
            if not symbolic or isinstance(I, sym.Integral):
                integrand = sym.lambdify([x], integrand)
                I = mpmath.quad(integrand, [Omega[0], Omega[1]])
            A[i,j] = A[j,i] = I

        integrand = psi[i]*f
        if symbolic:
            I = sym.integrate(integrand, (x, Omega[0], Omega[1]))
        if not symbolic or isinstance(I, sym.Integral):
            integrand = sym.lambdify([x], integrand)
            I = mpmath.quad(integrand, [Omega[0], Omega[1]])
        b[i,0] = I
    c = A.LUsolve(b)  # symbolic solve
    # c is a sympy Matrix object, numbers are in c[i,0]
    c = [sym.simplify(c[i,0]) for i in range(c.shape[0])]
    u = sum(c[i]*psi[i] for i in range(len(psi)))
    return u, c

def comparison_plot(f, u, Omega):
    x = sym.Symbol('x')
    f = sym.lambdify([x], f, modules="numpy")
    u = sym.lambdify([x], u, modules="numpy")
    resolution = 401  # no of points in plot
    xcoor  = np.linspace(Omega[0], Omega[1], resolution)
    exact  = f(xcoor)
    approx = u(xcoor)
    plt.plot(xcoor, approx)
    plt.plot(xcoor, exact)
    plt.legend(['approximation', 'exact'])
    plt.show()

x = sym.Symbol('x')
f = 10*(x-1)**2-1
Omega = [1, 2]
u, c = least_squares(f=f, psi=[1, x], Omega=Omega)

comparison_plot(f, u, Omega)

enter image description here

But if I replace the function with a complicated piecewise function: enter image description here

I got an error when I ran the same Python code again:

x = sym.Symbol('x')

f = abs(x)*sym.atan(x)/4 + x**2*sym.sin(1/x)
Omega = [-1, 1]
u, c = least_squares(f=f, psi=[1, x], Omega=Omega, symbolic=False)

comparison_plot(f, u, Omega)

Error is as follows:

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-68-04753e90fafc> in <module>
     51 f = abs(x)*sym.atan(x)/4 + x**2*sym.sin(1/x)
     52 Omega = [-1, 1]
---> 53 u, c = least_squares(f=f, psi=[1, x], Omega=Omega, symbolic=False)
     54 
     55 comparison_plot(f, u, Omega)

<ipython-input-68-04753e90fafc> in least_squares(f, psi, Omega, symbolic)
     26         if not symbolic or isinstance(I, sym.Integral):
     27             integrand = sym.lambdify([x], integrand)
---> 28             I = quad(integrand, [Omega[0], Omega[1]])
     29         b[i,0] = I
     30     c = A.LUsolve(b)  # symbolic solve

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in quad(ctx, f, *points, **kwargs)
    743             ctx.prec += 20
    744             if dim == 1:
--> 745                 v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
    746             elif dim == 2:
    747                 v, err = rule.summation(lambda x: \

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in summation(self, f, points, prec, epsilon, max_degree, verbose)
    231                     print("Integrating from %s to %s (degree %s of %s)" % \
    232                         (ctx.nstr(a), ctx.nstr(b), degree, max_degree))
--> 233                 result = self.sum_next(f, nodes, degree, prec, results, verbose)
    234                 results.append(result)
    235                 if degree > 1:

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in sum_next(self, f, nodes, degree, prec, previous, verbose)
    305         else:
    306             S = self.ctx.zero
--> 307         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    308         return h*S
    309 

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/ctx_mp_python.py in fdot(ctx, A, B, conjugate)
    934         hasattr_ = hasattr
    935         types = (ctx.mpf, ctx.mpc)
--> 936         for a, b in A:
    937             if type(a) not in types: a = ctx.convert(a)
    938             if type(b) not in types: b = ctx.convert(b)

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in <genexpr>(.0)
    305         else:
    306             S = self.ctx.zero
--> 307         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    308         return h*S
    309 

<lambdifygenerated-58> in _lambdifygenerated(x)
      1 def _lambdifygenerated(x):
----> 2     return x**2*sin(x**(-1.0)) + (1/4)*abs(x)*arctan(x)

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/ctx_mp_python.py in __pow__(self, other)

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/libmp/libelefun.py in mpf_pow(s, t, prec, rnd)
    326         raise ComplexResult("negative number raised to a fractional power")
    327     if texp >= 0:
--> 328         return mpf_pow_int(s, (-1)**tsign * (tman<<texp), prec, rnd)
    329     # s**(n/2) = sqrt(s)**n
    330     if texp == -1:

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/libmp/libmpf.py in mpf_pow_int(s, n, prec, rnd)
   1068         bc += bctable[int(man>>bc)]
   1069         return normalize1(0, man, exp+exp, bc, prec, rnd)
-> 1070     if n == -1: return mpf_div(fone, s, prec, rnd)
   1071     if n < 0:
   1072         inverse = mpf_pow_int(s, -n, prec+5, reciprocal_rnd[rnd])

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/libmp/libmpf.py in mpf_div(s, t, prec, rnd)
    958             return fzero
    959         if t == fzero:
--> 960             raise ZeroDivisionError
    961         s_special = (not sman) and sexp
    962         t_special = (not tman) and texp

ZeroDivisionError:

I don't know how to fix it.

Thank you for your suggestions, I wrote the function as a piecewise function as you said

from sympy import Piecewise, Symbol, sin, And, atan
f_1 = abs(x)*atan(x)/4 + x**2*sin(1/x)
f = Piecewise((f_1, And(-1 <= x) & (x < 0) | (x > 0) & (x <= 1)), (0, True))

The integral interval was also modified as follows：

I = quad(integrand, [Omega[0], 0, Omega[1]])

But something else went wrong:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
AttributeError: 'mpf' object has no attribute 'sin'

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
<ipython-input-149-344dd158c52d> in <module>
     54 
     55 Omega = [-1, 1]
---> 56 u, c = least_squares(f=f, psi=[1, x], Omega=Omega, symbolic=False)
     57 
     58 comparison_plot(f, u, Omega)

<ipython-input-149-344dd158c52d> in least_squares(f, psi, Omega, symbolic)
     27         if not symbolic or isinstance(I, sym.Integral):
     28             integrand = sym.lambdify([x], integrand)
---> 29             I = quad(integrand, [Omega[0], 0, Omega[1]])
     30         b[i,0] = I
     31     c = A.LUsolve(b)  # symbolic solve

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in quad(ctx, f, *points, **kwargs)
    743             ctx.prec += 20
    744             if dim == 1:
--> 745                 v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
    746             elif dim == 2:
    747                 v, err = rule.summation(lambda x: \

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in summation(self, f, points, prec, epsilon, max_degree, verbose)
    231                     print("Integrating from %s to %s (degree %s of %s)" % \
    232                         (ctx.nstr(a), ctx.nstr(b), degree, max_degree))
--> 233                 result = self.sum_next(f, nodes, degree, prec, results, verbose)
    234                 results.append(result)
    235                 if degree > 1:

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in sum_next(self, f, nodes, degree, prec, previous, verbose)
    305         else:
    306             S = self.ctx.zero
--> 307         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    308         return h*S
    309 

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/ctx_mp_python.py in fdot(ctx, A, B, conjugate)
    934         hasattr_ = hasattr
    935         types = (ctx.mpf, ctx.mpc)
--> 936         for a, b in A:
    937             if type(a) not in types: a = ctx.convert(a)
    938             if type(b) not in types: b = ctx.convert(b)

/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/mpmath/calculus/quadrature.py in <genexpr>(.0)
    305         else:
    306             S = self.ctx.zero
--> 307         S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
    308         return h*S
    309 

<lambdifygenerated-70> in _lambdifygenerated(x)
      1 def _lambdifygenerated(x):
----> 2     return select([logical_or.reduce((logical_and.reduce((greater_equal(x, -1),less(x, 0))),logical_and.reduce((less_equal(x, 1),greater(x, 0))))),True], [x**2*sin(x**(-1.0)) + (1/4)*abs(x)*arctan(x),0], default=nan)

TypeError: loop of ufunc does not support argument 0 of type mpf which has no callable sin method
Solution

  • According to this link: https://mpmath.org/doc/current/calculus/integration.html

    at "Highly variable functions":

    Whether splitting the interval or increasing the degree is more efficient differs from case to case. Another example is the function 1/(1+x2), which has a sharp peak centered around x=0:

    f = lambda x: 1/(1+x**2)
quad(f, [-100, 100])   # Bad
3.64804647105268
quad(f, [-100, 100], maxdegree=10)   # Good
3.12159332021646
quad(f, [-100, 0, 100])   # Also good
3.12159332021646