failure to verify linearity of l(v) = f*v when f is a polynomial

[campospinto]

Here is a small code which breaks (for f2 and degree = 5):

from sympde.expr.expr import LinearForm, integral
from sympde.topology                    import Square, ScalarFunctionSpace
from sympde.topology  import element_of


for degree in [3,4,5]:
        
    domain = Square()
    V = ScalarFunctionSpace('V', domain)
    v = element_of(V, name='v')

    x, y = domain.coordinates      
    f1 = (x-1)**degree * (y-1.1)**(2)*(y-1)**(degree-2)

    LinearForm(v, integral(domain, f1*v))
    print(f'f1, degree={degree}: OK')        

    f2 = (x-1.1)**degree * (y-1.1)**(2)*(y-1)**(degree-2)
    LinearForm(v, integral(domain, f2*v))
    print(f'f2, degree={degree}: OK')

[ratnania]

The problem is the numeric evaluation of (x-1.1)**degree where degree=5. You can see this in the following code:

def test_issue_141():
    from sympde.expr.expr import LinearForm
    from sympde.topology  import Square, ScalarFunctionSpace
    from sympde.topology  import element_of
    from sympde.expr.expr import is_linear_expression

    domain = Square()
    V = ScalarFunctionSpace('V', domain)
    v = element_of(V, name='v')

    x, y = domain.coordinates

    degree = 5
    g = lambda s,t: (x-s)**degree * (y-t)**(2)*(y-1)**(degree-2)
    f = lambda v,s,t : g(s,t) * v

    for s in [1, 1.1]:
        for t in [1, 1.1]:
            result = is_linear_expression(f(v,s,t), (v,), debug=False)
            print('>>> (s={s},t={t})'.format(s=s,t=t))
            print('    is linear = ', result)

The result is

>>> (s=1,t=1)
    is linear =  True
>>> (s=1,t=1.1)
    is linear =  True
>>> (s=1.1,t=1)
    is linear =  True
>>> (s=1.1,t=1.1)
    is linear =  False

There are two different ways to solve this issue:

• in your code: use a Constant instead of the actual value 1.1
• improve the code in is_linear_expression when testing the addition we do this:

    # ... check addition
    newargs = [left + right for left, right in zip(left_args, right_args)]

    newexpr    = expr.subs(zip(args, newargs))
    left_expr  = expr.subs(zip(args, left_args))
    right_expr = expr.subs(zip(args, right_args))

    a = newexpr
    b = left_expr + right_expr

    if not( (a-b).expand() == 0 or a.expand() == b.expand()):
        # TODO use a warning or exception?
        if debug:
            print('Failed to assert addition property')
            print('{} != {}'.format(a.expand(), b.expand()))
        return False
    # ...

we should check equality up to a given tolerance. But this needs further investigation since we're using symbolic expressions.