Why is SymPy's linsolve function not able to solve for the second to the last variable in this problem?

Question

(TLDR: SymPy's linsolve function is unable to solve the system of linear equations, generated from applying the finite difference method to an ODE BVP, when you pass the equations as a pure Python list, but is able to do so putting said equations list inside SymPy's Matrix function. This could be a bug that needs to be fixed especially considering that in the SymPy documentation, the example that they give you has them passing a list as an argument to linsolve. )

I have a boundary-value-problem ordinary differential equation that I intend to solve using the finite difference method. My ODE, in SymPy representation, is x*y(x).diff(x,2) + y(x).diff(x) + 500 = 0, with y(1)=600 and y(3.5)=25. Entire code is as follows:

import sympy as sp
from numpy import *
from matplotlib.pyplot import *
y = sp.Function('y')
ti = time.time()
steps = 10
ys = [ y(i) for i in range(steps+1) ]
ys[0], ys[-1] = 600, 25
xs = linspace(1, 3.5, steps+1)
dx = xs[1] - xs[0]
eqns = [ xs[i]*(ys[i-1] - 2*ys[i] + ys[i+1])/dx**2 + 
        ( ys[i+1] - ys[i-1] )/2/dx + 500
        for i in range(1, steps) ]
ys[1:-1] = sp.linsolve(eqns, ys[1:-1]).args[0]
scatter(xs, ys)
tf = time.time()
print(f'Time elapsed: {tf-ti}')

For ten steps, this works just fine. However, if I go higher than that immediately like 11 steps, SymPy is no longer able to solve the system of linear equations. Trying to plot the results throws a TypeError: can't convert expression to float. Examining the list of y values ys reveals that one of the variables, specifically the second to the last one, wasn't being solved for by sympy's linsolve, but instead the solutions to the other variables are being solved in terms of this unsolved variable. For example, if I have 50 steps, the second-to-the last variable is y(49) and this is featured on the solutions of the other unknowns and not solved for e.g. 2.02061855670103*y(49) - 26.1340206185567. In contrast, in another BVP ODE that I solved, y(x).diff(x,2) + y(x).diff(x) + y(x) - 1 with y(0)=1.5 and y(3)=2.5, it doesn't have any issue whether I want 10, 50, or 200 steps. It solves all the variables just fine, but this seems to be a peculiar exception as I encountered the aforementioned issue with many other ODEs.

SymPy's inconsistency here was quite frustrating. The only consolation is that before I ran into this problem, I have actually solved it a few times already with the varying number of steps that I wanted. I had enclosed the eqns variable inside sympy's Matrix function as in ys[1:-1] = sp.linsolve(sp.Matrix(eqns), ys[1:-1]).args[0] simply because it displayed better that way in the terminal. But for solving in a script file, I thought that wrapping it inside sp.Matrix is unnecessary and I naturally removed it to simplify things.

Answer 1

It is polite when formatting a question for SO (or anywhere else) to provide a complete code example without missing imports etc. Also you should distil the problem down to the minimum case and remove all of the unnecessary details. With that in mind a better way to demonstrate the issue is by actually figuring out what the arguments to linsolve are and presenting them directly e.g.:

from sympy import *

y = Function('y')

eqns = [
    -47.52*y(1) + 25.96*y(2) + 13436.0,
    25.96*y(1) - 56.32*y(2) + 30.36*y(3) + 500,
    30.36*y(2) - 65.12*y(3) + 34.76*y(4) + 500,
    34.76*y(3) - 73.92*y(4) + 39.16*y(5) + 500,
    39.16*y(4) - 82.72*y(5) + 43.56*y(6) + 500,
    43.56*y(5) - 91.52*y(6) + 47.96*y(7) + 500,
    47.96*y(6) - 100.32*y(7) + 52.36*y(8) + 500,
    52.36*y(7) - 109.12*y(8) + 56.76*y(9) + 500,
    56.76*y(8) - 117.92*y(9) + 61.16*y(10) + 500,
    61.16*y(9) - 126.72*y(10) + 2139.0
]

syms = [y(1), y(2), y(3), y(4), y(5), y(6), y(7), y(8), y(9), y(10)]

print(linsolve(eqns, syms))

Here you hoped to get a simple numerical solution for each of the unknowns but instead the result returned (from SymPy 1.8) is:

FiniteSet((5.88050359812056*y(10) - 5.77315239260531, 10.7643116711359*y(10) - 528.13328974178, 14.9403214726998*y(10) - 991.258097567359, 9.85496358613721e+15*y(10) - 1.00932650309452e+18, 7.35110502818395*y(10) - 312.312287998229, 5.84605452313345*y(10) - 217.293922525318, 4.47908204606922*y(10) - 141.418192750506, 3.22698120573309*y(10) - 81.4678489766327, 2.07194244604317*y(10) - 34.9738391105298, 1.0*y(10)))

The linsolve function will return a solution involving one or more unknowns if the system does not have a unique solution. Note also that there are some large numbers like 9.85496358613721e+15 which suggests that there might be numerical problems here.

Actually this is a bug in SymPy and it has already been fixed on the master branch: https://github.com/sympy/sympy/pull/21527

If you install SymPy from git then you can find the following as output instead:

FiniteSet((596.496767861074, 574.326903264955, 538.901024315178, 493.575084012669, 440.573815245681, 381.447789421181, 317.320815173574, 249.033388036155, 177.23053946471, 102.418085492911))

Also note that it is generally better to avoid using floats in SymPy as it is a library that is designed for exact symbolic computation. Solving a system of floating point equations like this can be done much more efficiently using NumPy or some other fixed precision floating point library that can use BLAS/LAPACK routines. To use rational arithmetic with sympy here you just need to change your linspace line to

xs = [sp.Rational(x) for x in linspace(1, 3.5, steps+1)]

which will then work fine with SymPy 1.8 and is in fact faster (at least if you have gmpy2 installed).