import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import quad
from scipy.interpolate import interp1d
from scipy.optimize import curve_fit, fsolve
import sympy as sp
sp.init_printing()
from IPython.display import display
def fη(κ):
η = 1
cη = 1
return (np.exp(-5.2*((((κ*η)**4) + cη**4)**0.25 - cη))
I = quad(fη, 0,1)[0]
File "<ipython-input-125-db6ae7012ebd>", line 6 I = quad(fη, 0,1)[0] ^ SyntaxError: invalid syntax
quad
statement below.Above, there is a missing ) at the end.
Count the number of left ( and right ) and make sure they match.
Delete any parentheses that are not needed (around the whole expression, and around the (κ*η)**4
return np.exp(-5.2*(((κ*η)**4 + cη**4)**0.25 - cη))
Add spaces to make the expression more clear
return np.exp( -5.2*( ((κ*η)**4 + cη**4)**0.25 - cη ) )
Compare carefully to the intended equation. $$f_{\eta}(\kappa) = \exp(-5.2([(\kappa\eta)^4 + c_{\eta}^4]^{1/4}-c_{\eta})).$$
Use sympy to solve for the root of $$x^3 + 15x^2 = 3x - 10$$ Then numerically evaluate each root.
x = sp.symbols('x')
ex = sp.Eq(x**3 + 15*x*x, 3*x - 10)
roots = sp.solve(ex,x)
display(roots)
roots.evalf()
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-126-f9b81b4cfd78> in <module> 3 roots = sp.solve(ex,x) 4 display(roots) ----> 5 roots.evalf() AttributeError: 'list' object has no attribute 'evalf'
roots
is a list and Python lists don't know anything about sympy's evalf
function.list.evalf()
.some_expression.evalf()
does make sense.Consider this problem and its solution.
Solve the following two equations in the two unknowns $c_L$ and $c_{\eta}$: $$\int_0^{\infty}E(\kappa)d\kappa = K,$$
$$\int_0^{\infty}2\nu\kappa^2E(\kappa)d\kappa = \epsilon.$$Here, $\nu = 1.5\times 10^{-5}$, $K = 0.008$, and $\epsilon = 0.003$.
You can take $\infty\approx 1E5$.
We also have the following functions and relations: $$L = \frac{K^{3/2}}{\epsilon},$$ $$\eta = \left(\frac{\nu^3}{\epsilon}\right)^{1/4},$$
$$E(\kappa) = 1.5\epsilon^{2/3}\kappa^{-5/3}f_L(\kappa)f_{\eta}(\kappa),$$$$f_L(\kappa) = \left(\frac{\kappa L}{[(\kappa L)^2 + c_L]^{1/2}}\right)^{11/3},$$$$f_{\eta}(\kappa) = \exp(-5.2([(\kappa\eta)^4 + c_{\eta}^4]^{1/4}-c_{\eta})).$$def F(cLcη):
cL = cLcη[0]
cη = cLcη[1]
ν = 1.5E-5
K = 0.008
ϵ = 0.003
inf = 1E5
L = K**1.5/ϵ
η = (ν**3/ϵ)**0.25
def fL(κ):
return (κ*L/np.sqrt((κ*L)**2 + cL))**(11/3)
def fη(κ):
return np.exp(-5.2*(((κ*η)**4 + cη**4)**0.25) - cη)
def E(κ):
return 1.5*ϵ**(2/3)*κ**(-5/3)*fL(κ)*fη(κ)
def twoNuK2E(κ):
return 2*ν*κ*2*E(κ)
eq0 = quad(E, 0, inf) - K
eq1 = quad(twoNuK2E, 0, inf)[0] - ϵ
return np.array([eq0, eq1])
#-----------------------------------------------------------
guess = np.array([5, 0.5])
cLcη = fsolve(F, guess)
print(f'cL = {cLcη[0]:.4f}')
print(f'cη = {cLcη[1]:.4f}')
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-127-871980a29428> in <module> 32 guess = np.array([5, 0.5]) 33 ---> 34 cLcη = fsolve(F, guess) 35 36 print(f'cL = {cLcη[0]:.4f}') /usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag) 146 'diag': diag} 147 --> 148 res = _root_hybr(func, x0, args, jac=fprime, **options) 149 if full_output: 150 x = res['x'] /usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options) 212 if not isinstance(args, tuple): 213 args = (args,) --> 214 shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,)) 215 if epsfcn is None: 216 epsfcn = finfo(dtype).eps /usr/local/lib/python3.7/site-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape) 25 def _check_func(checker, argname, thefunc, x0, args, numinputs, 26 output_shape=None): ---> 27 res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) 28 if (output_shape is not None) and (shape(res) != output_shape): 29 if (output_shape[0] != 1): <ipython-input-127-871980a29428> in F(cLcη) 23 return 2*ν*κ*2*E(κ) 24 ---> 25 eq0 = quad(E, 0, inf) - K 26 eq1 = quad(twoNuK2E, 0, inf)[0] - ϵ 27 TypeError: unsupported operand type(s) for -: 'tuple' and 'float'
return np.exp(-5.2*(((κ*η)**4 + cη**4)**0.25) - cη)
should be
return np.exp(-5.2*(((κ*η)**4 + cη**4)**0.25 - cη))
^
instead of **
6.02*10**23
instead of 6.02E23 (more of a computational sin than a bug).f
for both a function and a variable.I = quad(f(x), a, b)[0]
I = quad(f, a, b)[0]