I have following system of nonlinear equations:
This is very similar to Cholesky’s decomposition, but unfortunately it’s not it.
I have found very familiar solution: Solve a system of non-linear equations in Python (scipy.optimize.fsolve) but I don't know how to dynamically setup all equations of my system.
How to solve this system within Numpy, Scipy or any other packages in Python?
SymPy is the package that you're looking for. Here is a function that dynamically sets up your equations as described in your system:
import sympy as sp
def not_choleskys_decomposition(N):
# Initialisation of list of variables.
c = []
K = []
for n in range(N + 1):
c.append(sp.Symbol("c_{}".format(n), real=True))
K.append(sp.Symbol("K_{}".format(n), real=True))
# Setup your N+1 equations.
equations = []
for n in range(N + 1):
num_c = N + 1 - n
c_terms = [c_i * c_j for c_i, c_j in zip(c[:num_c], c[n:][:num_c])]
lhs = sp.Add(*c_terms)
equations.append(sp.Eq(lhs, K[n]))
return equations, c, K
Printing for verification can be done using the sympy.latex function which converts the math directly to LaTeX or the sympy.pprint function.
>>> test, _, _ = not_choleskys_decomposition(4)
>>> for t in test:
... sp.pprint(t)
c₀⋅c₀ + c₁⋅c₁ + c₂⋅c₂ + c₃⋅c₃ + c₄⋅c₄ = K₀
c₀⋅c₁ + c₁⋅c₂ + c₂⋅c₃ + c₃⋅c₄ = K₁
c₀⋅c₂ + c₁⋅c₃ + c₂⋅c₄ = K₂
c₀⋅c₃ + c₁⋅c₄ = K₃
c₀⋅c₄ = K₄
Replacing your variables can be done with the .subs method which belongs to each sympy.Eq object:
>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> dict_replace = {c_sym[0]: 1.0}
>>> new_equation = equations[0].subs(dict_replace)
>>> sympy.pprint(new_equation)
c₁⋅c₂ + 1.0⋅c₁ + c₂⋅c₃ + c₃⋅c₄ = K₁
Here's an example function that replaces your list of c or K.
def substitute_sym(equations, sym_list, val_list):
# Confirm all dimensions are consistent.
try:
assert len(equations) == len(sym_list) == len(val_list)
except AssertionError:
raise IndexError("Inconsistent dimensions.")
# Replace c symbols with values in equations.
substituted_equations = []
for eq in equations:
_eq = eq.subs(dict(zip(sym_list, val_list)))
substituted_equations.append(_eq)
return substituted_equations
Testing it on c:
>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> test_2 = substitute_sym(test_1, c_sym, [1] * 5)
>>> for t in test_2:
... sp.pprint(t)
5 = K₀
4 = K₁
3 = K₂
2 = K₃
1 = K₄
Testing it on K:
>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> K_vals = [1, 1.4, 0.2, 0.5, 3.0]
>>> test_3 = substitute_sym(test_1, K_sym, K_vals)
>>> for t in test_3:
... sp.pprint(t)
c₀⋅c₀ + c₁⋅c₁ + c₂⋅c₂ + c₃⋅c₃ + c₄⋅c₄ = 1
c₀⋅c₁ + c₁⋅c₂ + c₂⋅c₃ + c₃⋅c₄ = 1.4
c₀⋅c₂ + c₁⋅c₃ + c₂⋅c₄ = 0.2
c₀⋅c₃ + c₁⋅c₄ = 0.5
c₀⋅c₄ = 3.0
Solving your set of equation can be done with sympy.solvers.solveset.nonlinsolve(system, *symbols) documented here.
Solving for c given K
>>> from sympy.solvers.solveset import nonlinsolve
>>> test_solve, c_sym, K_sym = not_choleskys_decomposition(1)
>>> test_replaced = substitute_sym(test_solve, K_sym, [sp.Rational(0.5)] * 2)
>>> solved = nonlinsolve(test_replaced, c_sym)
>>> sp.pprint(solved)
⎧⎛ ⎛ 2⎞ ⎞ ⎛ ⎛
⎪⎜ ⎜ ⎛ √6 √2⋅ⅈ⎞ ⎟ ⎛ √6 √2⋅ⅈ⎞ √6 √2⋅ⅈ⎟ ⎜ ⎜ ⎛ √6 √2⋅ⅈ
⎨⎜-⎜-1 + 2⋅⎜- ── - ────⎟ ⎟⋅⎜- ── - ────⎟, - ── - ────⎟, ⎜-⎜-1 + 2⋅⎜- ── + ────
⎪⎝ ⎝ ⎝ 4 4 ⎠ ⎠ ⎝ 4 4 ⎠ 4 4 ⎠ ⎝ ⎝ ⎝ 4 4
⎩
2⎞ ⎞ ⎛ ⎛ 2⎞
⎞ ⎟ ⎛ √6 √2⋅ⅈ⎞ √6 √2⋅ⅈ⎟ ⎜ ⎜ ⎛√6 √2⋅ⅈ⎞ ⎟ ⎛√6 √2⋅ⅈ⎞ √6 √2⋅
⎟ ⎟⋅⎜- ── + ────⎟, - ── + ────⎟, ⎜-⎜-1 + 2⋅⎜── - ────⎟ ⎟⋅⎜── - ────⎟, ── - ───
⎠ ⎠ ⎝ 4 4 ⎠ 4 4 ⎠ ⎝ ⎝ ⎝4 4 ⎠ ⎠ ⎝4 4 ⎠ 4 4
⎞ ⎛ ⎛ 2⎞ ⎞⎫
ⅈ⎟ ⎜ ⎜ ⎛√6 √2⋅ⅈ⎞ ⎟ ⎛√6 √2⋅ⅈ⎞ √6 √2⋅ⅈ⎟⎪
─⎟, ⎜-⎜-1 + 2⋅⎜── + ────⎟ ⎟⋅⎜── + ────⎟, ── + ────⎟⎬
⎠ ⎝ ⎝ ⎝4 4 ⎠ ⎠ ⎝4 4 ⎠ 4 4 ⎠⎪
⎭
I hope this helps get you jump started with SymPy if you choose to use it. I have rarely used its nonlinear solver so I cannot guarantee anything past the examples. Arguments are also relevant and solvers are temperamental. You'll notice I replaced K with sp.Rational(0.5) for each index. It was done to avoid an error that is thrown otherwise when defining float type for some solvers. Good luck.
EDIT:
Also note that you don't need to use the solvers in SymPy. I rarely do. I use the package for the symbolic mathematics in Python for LaTeX and equation manipulation.
EDIT:
Making it work with scipy.optimize.fsolve
from sympy.utilities.lambdify import lambdify
from scipy.optimize import fsolve
def lambdify_equations_for_solving_c(equations, c_sym, K_sym, K_val):
f = []
equations_subbed = substitute_sym(equations, K_sym, K_val)
for eq in equations_subbed:
# Note that I reformulate equation into an expression here
# for determining the roots.
f.append(lambdify(c_sym, eq.args[0] - eq.args[1]))
return f
N = 4
equations, c_sym, K_sym = not_choleskys_decomposition(N)
K_vals = [0.1, 0, 0.1, 0, 0]
f = lambdify_equations_for_solving_c(equations, c_sym, K_sym, K_vals)
def fsolve_friendly(p):
return [_f(*p) for _f in f]
c_sol = fsolve(fsolve_friendly, x0=(0, 0.1, -1.1, 0, 0))
fsolve returns warnings for non-converged results.
/home/ggarrett/anaconda3/envs/sigh/lib/python3.7/site-packages/scipy/optimize/minpack.py:162: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last five Jacobian evaluations.
warnings.warn(msg, RuntimeWarning)
That's a complete answer on how to form your equations to be solved as python callables. The solver routine is another challenge itself.
