1

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
Improve this question
8
  • @DavisHerring Thanks for your reminder, it has been revised. Commented Oct 8, 2022 at 7:12
  • 1
    It looks like the domain in your code includes zero: Omega = [-1, 1] Commented Oct 8, 2022 at 7:15
  • I'm not familiar with sympy, but it looks like there is something called Piecewise; maybe you can define a piecewise domain with that: stackoverflow.com/questions/39030227/… Commented Oct 8, 2022 at 7:17
  • As anon01 says, your domain includes 0. There is a 1/x part in your function which is resulting a zero division error. We can revise: the lambda function, the range or the syntax of quad. Commented Oct 8, 2022 at 7:30
  • @anon01 I wrote the function as a piecewise function as you said, but something else went wrong, I wrote it in my problem. Commented Oct 8, 2022 at 8:37

1 Answer 1

Reset to default
1

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
Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

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))
but now you have some other error: AttributeError Traceback (most recent call last) AttributeError: 'mpf' object has no attribute 'sin'
man.hubwiz.com/docset/SymPy.docset/Contents/Resources/Documents/… | lambdify can receive a third parameter | search and look around this line >>> f = lambdify(x, sin(x), [{'sin': mysin}, 'numpy']) | below is another one >>> g = lambdify(x, x + sin(x), 'numpy')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.