I am attempting to solve an equation utilising Sympy, in combination with Uncertainties to perform error propagation. However, due to the nature of the equation, the exact form of the solution is a complex number which causes problems when trying to acquire a numerical solution.
Here's a MWE of my code:
from uncertainties import ufloat
import uncertainties.umath as um
import sympy as sp
from sympy.utilities.lambdify import lambdify
inc, ratio = sp.symbols('inc ratio')
beta = sp.Symbol("beta")
eqn=sp.Eq(sp.cos(inc+beta)/sp.cos(inc-beta),ratio)
s = sp.solve(eqn,beta)[1] # unsure if [0] or [1] is the best solution?
sol = lambdify([inc, ratio],s,modules=um) # necessary for integrating ufloat inputs as an option
i = ufloat(0.923, 0.1202)
r = ufloat(4.089, 0.1844)
b = sol(i,r) # here is the problem!!!
This returns
TypeError: unsupported operand type(s) for *: 'complex' and 'AffineScalarFunc'
My outputs for s and help(sol) are:
[-I*log(-(ratio*exp(2*I*inc) - 1)/(ratio - exp(2*I*inc)))/2,
-I*log(-sqrt(-(ratio*exp(2*I*inc) - 1)/(ratio - exp(2*I*inc))))]
Help on function _lambdifygenerated in module uncertainties.umath:
_lambdifygenerated(inc, ratio)
Created with lambdify. Signature:
func(inc, ratio)
Expression:
[-I*log(-(ratio*exp(2*I*inc) - 1)/(ratio - exp(2*I*inc)))/2,...
Source code:
def _lambdifygenerated(inc, ratio):
return [-1/2*1j*log(-(ratio*exp(2*1j*inc) - 1)/(ratio - exp(2*1j*inc))), -1j*log(-sqrt(-(ratio*exp(2*1j*inc) - 1)/(ratio - exp(2*1j*inc))))]
Imported modules:
This question has been very helpful with getting the two packages to cooperate together, but unfortunately doesn't help with solving my particular problem. I know the equation has a real solution (likely something modulo 2pi). I have tried constraining my variables to be only real, and that doesn't help either. Thanks in advance!
Could you use the symbolic propogation of errors formula to find the result?
# your code then...
>>> from sympy import re, sqrt
>>> usym = {v: Symbol('_'+v.name) for v in s.free_symbols}
>>> u2=sum([j**2 for j in [s.diff(v)*usym[v] for v in s.free_symbols]])
>>> reps = {inc: 0.923, _inc:.1202, ratio:4.089, _ratio:.1844}
>>> snum = s.subs(reps).n()
>>> unum = sqrt(re(u2.subs(reps))).n()
>>> ufloat(snum, unum)
2.711004137405918+/-0.09516734706025402