CIE_4140_Lecture_3_4_Python#
Integral Fourier Transform and frequency-domain analysis#
import sympy as sp
tau= sp.symbols('tau',real=True,positive=True)
t = sp.symbols('t',real=True)
Example Function of Time
f = sp.Heaviside(t-tau)-sp.Heaviside(t-2*tau)
sp.plot(f.subs(tau,2),(t,0,10));
omega = sp.symbols('omega',real=True,positive=True)
phi = sp.symbols('phi',real=True,positive=True)
f_omega = sp.integrate(f*sp.exp(-sp.I*omega*t), (t,-sp.oo,sp.oo))
f_phi = sp.fourier_transform(f,t,phi)
f_omega = sp.simplify(f_omega)
f_phi = sp.simplify(f_phi)
display(f_omega)
display(f_phi)
\[\displaystyle \frac{i \left(1 - e^{i \omega \tau}\right) e^{- 2 i \omega \tau}}{\omega}\]
\[\displaystyle \frac{i \left(1 - e^{2 i \pi \phi \tau}\right) e^{- 4 i \pi \phi \tau}}{2 \pi \phi}\]
Another definition for the fourier transformation is used by sympy
p0 = sp.plotting.plot(sp.re(f_omega.subs([(tau,2)])),(omega,-10,10),label='Real part' ,legend=True,show=False,adaptive=False,nb_of_points=3000)
p1 = sp.plotting.plot(sp.im(f_omega.subs([(tau,2)])),(omega,-10,10),label='Imaginary part',legend=True,show=False,adaptive=False,nb_of_points=3000)
p0.append(p1[0])
p0.show()
p0 = sp.plotting.plot(sp.re(f_phi.subs([(tau,2)])),(phi,-10/2/sp.pi,10/2/sp.pi),label='Real part' ,legend=True,show=False,adaptive=False,nb_of_points=3000)
p1 = sp.plotting.plot(sp.im(f_phi.subs([(tau,2)])),(phi,-10/2/sp.pi,10/2/sp.pi),label='Imaginary part',legend=True,show=False,adaptive=False,nb_of_points=3000)
p0.append(p1[0])
p0.show()
x = sp.symbols('x',cls = sp.Function)
m, c, k = sp.symbols('m, c, k',real=True)#,positive=True)
Equation_of_Motion = m*sp.diff(x(t),t,2)+c*sp.diff(x(t),t)+k*x(t)-f
display(Equation_of_Motion)
\[\displaystyle c \frac{d}{d t} x{\left(t \right)} + k x{\left(t \right)} + m \frac{d^{2}}{d t^{2}} x{\left(t \right)} + \theta\left(t - 2 \tau\right) - \theta\left(t - \tau\right)\]
sp.fourier_transform(Equation_of_Motion,t,phi)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_21888\1221906460.py in <module>
----> 1 sp.fourier_transform(Equation_of_Motion,t,phi)
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in fourier_transform(f, x, k, **hints)
2258 mellin_transform, laplace_transform
2259 """
-> 2260 return FourierTransform(f, x, k).doit(**hints)
2261
2262
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in doit(self, **hints)
176 if fn.is_Add:
177 hints['needeval'] = needeval
--> 178 res = [self.__class__(*([x] + list(self.args[1:]))).doit(**hints)
179 for x in fn.args]
180 extra = []
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in <listcomp>(.0)
176 if fn.is_Add:
177 hints['needeval'] = needeval
--> 178 res = [self.__class__(*([x] + list(self.args[1:]))).doit(**hints)
179 for x in fn.args]
180 extra = []
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in doit(self, **hints)
169 hints['simplify'] = simplify
170
--> 171 fn, T = self._try_directly(**hints)
172
173 if T is not None:
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in _try_directly(self, **hints)
133 if try_directly:
134 try:
--> 135 T = self._compute_transform(self.function,
136 self.function_variable, self.transform_variable, **hints)
137 except IntegralTransformError:
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in _compute_transform(self, f, x, k, **hints)
2190
2191 def _compute_transform(self, f, x, k, **hints):
-> 2192 return _fourier_transform(f, x, k,
2193 self.a(), self.b(),
2194 self.__class__._name, **hints)
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in wrapper(noconds, *args, **kwargs)
253 @wraps(func)
254 def wrapper(*args, noconds=default, **kwargs):
--> 255 res = func(*args, **kwargs)
256 if noconds:
257 return res[0]
~\Anaconda3\lib\site-packages\sympy\integrals\transforms.py in _fourier_transform(f, x, k, a, b, name, simplify)
2164 return _simplify(F, simplify), S.true
2165
-> 2166 integral_f = integrate(f, (x, S.NegativeInfinity, S.Infinity))
2167 if integral_f in (S.NegativeInfinity, S.Infinity, S.NaN) or integral_f.has(Integral):
2168 raise IntegralTransformError(name, f, 'function not integrable on real axis')
~\Anaconda3\lib\site-packages\sympy\integrals\integrals.py in integrate(meijerg, conds, risch, heurisch, manual, *args, **kwargs)
1565
1566 if isinstance(integral, Integral):
-> 1567 return integral.doit(**doit_flags)
1568 else:
1569 new_args = [a.doit(**doit_flags) if isinstance(a, Integral) else a
~\Anaconda3\lib\site-packages\sympy\integrals\integrals.py in doit(self, **hints)
709 try:
710 evalued = Add(*others)._eval_interval(x, a, b)
--> 711 evalued_pw = piecewise_fold(Add(*piecewises))._eval_interval(x, a, b)
712 function = uneval + evalued + evalued_pw
713 except NotImplementedError:
~\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in _eval_interval(self, sym, a, b, _first)
648 # TODO simplify hi <= upto
649 return Piecewise((sum, hi <= upto), (Undefined, True))
--> 650 sum += abei[i][-2]._eval_interval(x, a, b)
651 upto = b
652 return sum
~\Anaconda3\lib\site-packages\sympy\core\expr.py in _eval_interval(self, x, a, b)
936 return S.Zero
937
--> 938 A = _eval_endpoint(left=True)
939 if A is S.NaN:
940 return A
~\Anaconda3\lib\site-packages\sympy\core\expr.py in _eval_endpoint(left)
925 S.ComplexInfinity, AccumBounds):
926 if (a < b) != False:
--> 927 C = limit(self, x, c, "+" if left else "-")
928 else:
929 C = limit(self, x, c, "-" if left else "+")
~\Anaconda3\lib\site-packages\sympy\series\limits.py in limit(e, z, z0, dir)
62 """
63
---> 64 return Limit(e, z, z0, dir).doit(deep=False)
65
66
~\Anaconda3\lib\site-packages\sympy\series\limits.py in doit(self, **hints)
307 newe = e.subs(z, z + z0)
308 try:
--> 309 coeff, ex = newe.leadterm(z, cdir=cdir)
310 except (ValueError, NotImplementedError, PoleError):
311 # The NotImplementedError catching is for custom functions
~\Anaconda3\lib\site-packages\sympy\core\expr.py in leadterm(self, x, logx, cdir)
3505 from .symbol import Dummy
3506 from sympy.functions.elementary.exponential import log
-> 3507 l = self.as_leading_term(x, logx=logx, cdir=cdir)
3508 d = Dummy('logx')
3509 if l.has(log(x)):
~\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
68 def wrapper(*args, **kwargs):
69 try:
---> 70 retval = cfunc(*args, **kwargs)
71 except TypeError as e:
72 if not e.args or not e.args[0].startswith('unhashable type:'):
~\Anaconda3\lib\site-packages\sympy\core\expr.py in as_leading_term(self, logx, cdir, *symbols)
3468 if x not in self.free_symbols:
3469 return self
-> 3470 obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
3471 if obj is not None:
3472 from sympy.simplify.powsimp import powsimp
~\Anaconda3\lib\site-packages\sympy\core\mul.py in _eval_as_leading_term(self, x, logx, cdir)
1971
1972 def _eval_as_leading_term(self, x, logx=None, cdir=0):
-> 1973 return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
1974
1975 def _eval_conjugate(self):
~\Anaconda3\lib\site-packages\sympy\core\mul.py in <listcomp>(.0)
1971
1972 def _eval_as_leading_term(self, x, logx=None, cdir=0):
-> 1973 return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
1974
1975 def _eval_conjugate(self):
~\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
68 def wrapper(*args, **kwargs):
69 try:
---> 70 retval = cfunc(*args, **kwargs)
71 except TypeError as e:
72 if not e.args or not e.args[0].startswith('unhashable type:'):
~\Anaconda3\lib\site-packages\sympy\core\expr.py in as_leading_term(self, logx, cdir, *symbols)
3468 if x not in self.free_symbols:
3469 return self
-> 3470 obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
3471 if obj is not None:
3472 from sympy.simplify.powsimp import powsimp
~\Anaconda3\lib\site-packages\sympy\functions\special\hyper.py in _eval_as_leading_term(self, x, logx, cdir)
692 def _eval_as_leading_term(self, x, logx=None, cdir=0):
693 from sympy.simplify.hyperexpand import hyperexpand
--> 694 return hyperexpand(self).as_leading_term(x, logx=logx, cdir=cdir)
695
696 def integrand(self, s):
~\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
68 def wrapper(*args, **kwargs):
69 try:
---> 70 retval = cfunc(*args, **kwargs)
71 except TypeError as e:
72 if not e.args or not e.args[0].startswith('unhashable type:'):
~\Anaconda3\lib\site-packages\sympy\core\expr.py in as_leading_term(self, logx, cdir, *symbols)
3468 if x not in self.free_symbols:
3469 return self
-> 3470 obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
3471 if obj is not None:
3472 from sympy.simplify.powsimp import powsimp
~\Anaconda3\lib\site-packages\sympy\functions\elementary\piecewise.py in _eval_as_leading_term(self, x, logx, cdir)
314 def _eval_as_leading_term(self, x, logx=None, cdir=0):
315 for e, c in self.args:
--> 316 if c == True or c.subs(x, 0) == True:
317 return e.as_leading_term(x)
318
~\Anaconda3\lib\site-packages\sympy\core\basic.py in subs(self, *args, **kwargs)
995 rv = self
996 for old, new in sequence:
--> 997 rv = rv._subs(old, new, **kwargs)
998 if not isinstance(rv, Basic):
999 break
~\Anaconda3\lib\site-packages\sympy\core\cache.py in wrapper(*args, **kwargs)
68 def wrapper(*args, **kwargs):
69 try:
---> 70 retval = cfunc(*args, **kwargs)
71 except TypeError as e:
72 if not e.args or not e.args[0].startswith('unhashable type:'):
~\Anaconda3\lib\site-packages\sympy\core\basic.py in _subs(self, old, new, **hints)
1109 rv = self._eval_subs(old, new)
1110 if rv is None:
-> 1111 rv = fallback(self, old, new)
1112 return rv
1113
~\Anaconda3\lib\site-packages\sympy\core\basic.py in fallback(self, old, new)
1086 args[i] = arg
1087 if hit:
-> 1088 rv = self.func(*args)
1089 hack2 = hints.get('hack2', False)
1090 if hack2 and self.is_Mul and not rv.is_Mul: # 2-arg hack
~\Anaconda3\lib\site-packages\sympy\core\relational.py in __new__(cls, lhs, rhs, **options)
832 for me in (lhs, rhs):
833 if me.is_extended_real is False:
--> 834 raise TypeError("Invalid comparison of non-real %s" % me)
835 if me is S.NaN:
836 raise TypeError("Invalid NaN comparison")
TypeError: Invalid comparison of non-real zoo
So inhomogeneous equation of motion cannot be Fourier transformed, homogeneous is possible:
Equation_of_Motion_in_frequency_domain_homogeneous = sp.fourier_transform(m*sp.diff(x(t),t,2)+c*sp.diff(x(t),t)+k*x(t),t,phi)
Equation_of_Motion_in_frequency_domain = Equation_of_Motion_in_frequency_domain_homogeneous - f_phi
display(Equation_of_Motion_in_frequency_domain)
\[\displaystyle c \mathcal{F}_{t}\left[\frac{d}{d t} x{\left(t \right)}\right]\left(\phi\right) + k \mathcal{F}_{t}\left[x{\left(t \right)}\right]\left(\phi\right) + m \mathcal{F}_{t}\left[\frac{d^{2}}{d t^{2}} x{\left(t \right)}\right]\left(\phi\right) - \frac{i \left(1 - e^{2 i \pi \phi \tau}\right) e^{- 4 i \pi \phi \tau}}{2 \pi \phi}\]
Derivatives of \({d \over {dt}}x\left( t \right)\) not evaluated
Equation_of_Motion_in_frequency_domain = Equation_of_Motion_in_frequency_domain.subs(sp.fourier_transform(x(t).diff(t),t,phi),sp.I*2*sp.pi*phi*sp.fourier_transform(x(t),t,phi))
Equation_of_Motion_in_frequency_domain = Equation_of_Motion_in_frequency_domain.subs(sp.fourier_transform(x(t).diff(t,2),t,phi),(sp.I*2*sp.pi*phi)**2*sp.fourier_transform(x(t),t,phi))
display(Equation_of_Motion_in_frequency_domain)
\[\displaystyle 2 i \pi c \phi \mathcal{F}_{t}\left[x{\left(t \right)}\right]\left(\phi\right) + k \mathcal{F}_{t}\left[x{\left(t \right)}\right]\left(\phi\right) - 4 \pi^{2} m \phi^{2} \mathcal{F}_{t}\left[x{\left(t \right)}\right]\left(\phi\right) - \frac{i \left(1 - e^{2 i \pi \phi \tau}\right) e^{- 4 i \pi \phi \tau}}{2 \pi \phi}\]
solution_in_frequency_domain = sp.solve(sp.Eq(Equation_of_Motion_in_frequency_domain,0),sp.FourierTransform(x(t), t, phi))[0]
display(solution_in_frequency_domain)
\[\displaystyle \frac{i \left(1 - e^{2 i \pi \phi \tau}\right) e^{- 4 i \pi \phi \tau}}{2 \pi \phi \left(2 i \pi c \phi + k - 4 \pi^{2} m \phi^{2}\right)}\]
p0 = sp.plotting.plot(sp.re(solution_in_frequency_domain.subs([(tau,2),(k,4),(m,1),(c,0.75)])),(phi,0.001,15/2/sp.pi),label='Real part' ,legend=True,show=False,adaptive=False,nb_of_points=3000)
p1 = sp.plotting.plot(sp.im(solution_in_frequency_domain.subs([(tau,2),(k,4),(m,1),(c,0.75)])),(phi,0.001,15/2/sp.pi),label='Imaginary part',legend=True,show=False,adaptive=False,nb_of_points=3000)
p0.append(p1[0])
p0.show()
solution = sp.inverse_fourier_transform(solution_in_frequency_domain, phi,t)
display(solution)
\[\displaystyle i \left(- \mathcal{F}^{-1}_{\phi}\left[\frac{1}{4 i \pi^{2} c \phi^{2} e^{2 i \pi \phi \tau} + 2 \pi k \phi e^{2 i \pi \phi \tau} - 8 \pi^{3} m \phi^{3} e^{2 i \pi \phi \tau}}\right]\left(t\right) + \mathcal{F}^{-1}_{\phi}\left[\frac{1}{4 i \pi^{2} c \phi^{2} e^{4 i \pi \phi \tau} + 2 \pi k \phi e^{4 i \pi \phi \tau} - 8 \pi^{3} m \phi^{3} e^{4 i \pi \phi \tau}}\right]\left(t\right)\right)\]
Above cell takes a long time, but is not fully solved.
solution.subs([(tau,2),(k,4),(m,1),(c,3/4),(t,2)]).evalf()
\[\displaystyle i \left(- \mathcal{F}^{-1}_{\phi}\left[\frac{1}{- 8 \pi^{3} \phi^{3} e^{4 i \pi \phi} + 3.0 i \pi^{2} \phi^{2} e^{4 i \pi \phi} + 8 \pi \phi e^{4 i \pi \phi}}\right]\left(2\right) + \mathcal{F}^{-1}_{\phi}\left[\frac{1}{- 8 \pi^{3} \phi^{3} e^{8 i \pi \phi} + 3.0 i \pi^{2} \phi^{2} e^{8 i \pi \phi} + 8 \pi \phi e^{8 i \pi \phi}}\right]\left(2\right)\right)\]
Function cannot be evaluated
Manual fourier transform is tried:
solution_1 = 2 * sp.integrate(sp.re(solution_in_frequency_domain*sp.exp(sp.I*2*sp.pi*phi*t)),(phi,0,sp.oo))
display(solution_1)
\[\displaystyle \frac{\int\limits_{0}^{\infty} \left(- \frac{\operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\right)\, d\phi + \int\limits_{0}^{\infty} \frac{\sin{\left(2 \pi \phi \tau \right)} \operatorname{re}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\, d\phi + \int\limits_{0}^{\infty} \frac{\cos{\left(2 \pi \phi \tau \right)} \operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\, d\phi}{\pi}\]
solution_1.subs([(tau,2),(k,4),(m,1),(c,3/4),(t,4)]).evalf()
\[\displaystyle 0.3493\]
Numerical evaluation seems possible, but takes too long. In Maple this takes quite a long time to solve as well.
solution_func = sp.lambdify(t,solution_1.subs([(tau,2),(k,4),(m,1),(c,3/4)]))
display(solution_func(4))
<lambdifygenerated-1>:2: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
return (quad(lambda phi: -imag(exp(2*1j*pi*phi*t)/(-4*pi**2*phi**2*exp(8*1j*pi*phi) + 1.5*1j*pi*phi*exp(8*1j*pi*phi) + 4*exp(8*1j*pi*phi)))/phi, 0, PINF)[0] + quad(lambda phi: sin(4*pi*phi)*real(exp(2*1j*pi*phi*t)/(-4*pi**2*phi**2*exp(8*1j*pi*phi) + 1.5*1j*pi*phi*exp(8*1j*pi*phi) + 4*exp(8*1j*pi*phi)))/phi, 0, PINF)[0] + quad(lambda phi: cos(4*pi*phi)*imag(exp(2*1j*pi*phi*t)/(-4*pi**2*phi**2*exp(8*1j*pi*phi) + 1.5*1j*pi*phi*exp(8*1j*pi*phi) + 4*exp(8*1j*pi*phi)))/phi, 0, PINF)[0])/pi
0.3493028443000069
import numpy as np
import matplotlib.pylab as plt
func_values = []
t_list = np.linspace(0,15,100)
for t_i in t_list:
func_values.append(solution_func(t_i))
plt.plot(t_list,func_values);
Another try with an upper limit of 10 instead of inifity allows evaluating of the function.
solution_2 = 2 * sp.integrate(sp.re(solution_in_frequency_domain*sp.exp(sp.I*2*sp.pi*phi*t)),(phi,0,1.5))
display(solution_2)
\[\displaystyle \frac{\int\limits_{0}^{1.5} \left(- \frac{\operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\right)\, d\phi + \int\limits_{0}^{1.5} \frac{\sin{\left(2 \pi \phi \tau \right)} \operatorname{re}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\, d\phi + \int\limits_{0}^{1.5} \frac{\cos{\left(2 \pi \phi \tau \right)} \operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{2 i \pi c \phi e^{4 i \pi \phi \tau} + k e^{4 i \pi \phi \tau} - 4 \pi^{2} m \phi^{2} e^{4 i \pi \phi \tau}}\right)}}{\phi}\, d\phi}{\pi}\]
solution_2.subs([(tau,2),(k,4),(m,1),(c,3/4),(t,4)]).evalf()
\[\displaystyle 0.349397566986233\]
sp.plot(solution_2.subs([(tau,2),(k,4),(m,1),(c,3/4)]),(t,0,15));
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_21888\313203652.py in <module>
----> 1 sp.plot(solution_2.subs([(tau,2),(k,4),(m,1),(c,3/4)]),(t,0,15));
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in plot(show, *args, **kwargs)
1871 plots = Plot(*series, **kwargs)
1872 if show:
-> 1873 plots.show()
1874 return plots
1875
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in show(self)
249 self._backend.close()
250 self._backend = self.backend(self)
--> 251 self._backend.show()
252
253 def save(self, path):
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in show(self)
1547
1548 def show(self):
-> 1549 self.process_series()
1550 #TODO after fixing https://github.com/ipython/ipython/issues/1255
1551 # you can uncomment the next line and remove the pyplot.show() call
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in process_series(self)
1544 if isinstance(self.parent, PlotGrid):
1545 parent = self.parent.args[i]
-> 1546 self._process_series(series, ax, parent)
1547
1548 def show(self):
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in _process_series(self, series, ax, parent)
1365 # Create the collections
1366 if s.is_2Dline:
-> 1367 x, y = s.get_data()
1368 if (isinstance(s.line_color, (int, float)) or
1369 callable(s.line_color)):
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in get_data(self)
603 """
604 np = import_module('numpy')
--> 605 points = self.get_points()
606 if self.steps is True:
607 if len(points) == 2:
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in get_points(self)
781 x_coords.append(self.start)
782 y_coords.append(f_start)
--> 783 sample(np.array([self.start, f_start]),
784 np.array([self.end, f_end]), 0)
785
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in sample(p, q, depth)
750 elif depth < 6:
751 sample(p, new_point, depth + 1)
--> 752 sample(new_point, q, depth + 1)
753
754 # Sample ten points if complex values are encountered
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in sample(p, q, depth)
749 # depth of 6. We are not using linspace to avoid aliasing.
750 elif depth < 6:
--> 751 sample(p, new_point, depth + 1)
752 sample(new_point, q, depth + 1)
753
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in sample(p, q, depth)
749 # depth of 6. We are not using linspace to avoid aliasing.
750 elif depth < 6:
--> 751 sample(p, new_point, depth + 1)
752 sample(new_point, q, depth + 1)
753
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in sample(p, q, depth)
750 elif depth < 6:
751 sample(p, new_point, depth + 1)
--> 752 sample(new_point, q, depth + 1)
753
754 # Sample ten points if complex values are encountered
~\Anaconda3\lib\site-packages\sympy\plotting\plot.py in sample(p, q, depth)
738 else:
739 xnew = p[0] + random * (q[0] - p[0])
--> 740 ynew = f(xnew)
741 new_point = np.array([xnew, ynew])
742
~\Anaconda3\lib\site-packages\sympy\plotting\experimental_lambdify.py in __call__(self, args)
174 try:
175 #The result can be sympy.Float. Hence wrap it with complex type.
--> 176 result = complex(self.lambda_func(args))
177 if abs(result.imag) > 1e-7 * abs(result):
178 return None
~\Anaconda3\lib\site-packages\sympy\plotting\experimental_lambdify.py in __call__(self, *args, **kwargs)
270
271 def __call__(self, *args, **kwargs):
--> 272 return self.lambda_func(*args, **kwargs)
273
274
<string> in <lambda>(x0)
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(self, n, subs, maxn, chop, strict, quad, verbose)
1657 options['quad'] = quad
1658 try:
-> 1659 result = evalf(self, prec + 4, options)
1660 except NotImplementedError:
1661 # Fall back to the ordinary evalf
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_integral(expr, prec, options)
1185 maxprec = options.get('maxprec', INF)
1186 while 1:
-> 1187 result = do_integral(expr, workprec, options)
1188 accuracy = complex_accuracy(result)
1189 if accuracy >= prec: # achieved desired precision
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in do_integral(expr, prec, options)
1146 quadrature_error = MINUS_INF
1147 else:
-> 1148 result, quadrature_err = quadts(f, [xlow, xhigh], error=1)
1149 quadrature_error = fastlog(quadrature_err._mpf_)
1150
~\Anaconda3\lib\site-packages\mpmath\calculus\quadrature.py in quadts(ctx, *args, **kwargs)
786 """
787 kwargs['method'] = 'tanh-sinh'
--> 788 return ctx.quad(*args, **kwargs)
789
790 def quadgl(ctx, *args, **kwargs):
~\Anaconda3\lib\site-packages\mpmath\calculus\quadrature.py in quad(ctx, f, *points, **kwargs)
743 ctx.prec += 20
744 if dim == 1:
--> 745 v, err = rule.summation(f, points[0], prec, epsilon, m, verbose)
746 elif dim == 2:
747 v, err = rule.summation(lambda x: \
~\Anaconda3\lib\site-packages\mpmath\calculus\quadrature.py in summation(self, f, points, prec, epsilon, max_degree, verbose)
231 print("Integrating from %s to %s (degree %s of %s)" % \
232 (ctx.nstr(a), ctx.nstr(b), degree, max_degree))
--> 233 result = self.sum_next(f, nodes, degree, prec, results, verbose)
234 results.append(result)
235 if degree > 1:
~\Anaconda3\lib\site-packages\mpmath\calculus\quadrature.py in sum_next(self, f, nodes, degree, prec, previous, verbose)
305 else:
306 S = self.ctx.zero
--> 307 S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
308 return h*S
309
~\Anaconda3\lib\site-packages\mpmath\ctx_mp_python.py in fdot(ctx, A, B, conjugate)
934 hasattr_ = hasattr
935 types = (ctx.mpf, ctx.mpc)
--> 936 for a, b in A:
937 if type(a) not in types: a = ctx.convert(a)
938 if type(b) not in types: b = ctx.convert(b)
~\Anaconda3\lib\site-packages\mpmath\calculus\quadrature.py in <genexpr>(.0)
305 else:
306 S = self.ctx.zero
--> 307 S += self.ctx.fdot((w,f(x)) for (x,w) in nodes)
308 return h*S
309
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in f(t)
1119 def f(t: 'Expr') -> tUnion[mpc, mpf]:
1120 nonlocal max_real_term, max_imag_term
-> 1121 re, im, re_acc, im_acc = evalf(func, mp.prec, {'subs': {x: t}})
1122
1123 have_part[0] = re or have_part[0]
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_mul(v, prec, options)
696 elif i == last and arg is S.One:
697 continue
--> 698 re, im, re_acc, im_acc = evalf(arg, working_prec, options)
699 if re and im:
700 complex_factors.append((re, im, re_acc, im_acc))
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_im(expr, prec, options)
311
312 def evalf_im(expr: 'im', prec: int, options: OPT_DICT) -> TMP_RES:
--> 313 return get_complex_part(expr.args[0], 1, prec, options)
314
315
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in get_complex_part(expr, no, prec, options)
291 i = 0
292 while 1:
--> 293 res = evalf(expr, workprec, options)
294 if res is S.ComplexInfinity:
295 return fnan, None, prec, None
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_mul(v, prec, options)
648 from .numbers import Float
649 for arg in args:
--> 650 result = evalf(arg, prec, options)
651 if result is S.ComplexInfinity:
652 special.append(result)
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_pow(v, prec, options)
776 # base must be evaluated with increased precision if p is large
777 prec += int(math.log(abs(p), 2))
--> 778 result = evalf(base, prec + 5, options)
779 if result is S.ComplexInfinity:
780 if p < 0:
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_add(v, prec, options)
600 options['maxprec'] = min(oldmaxprec, 2*prec)
601
--> 602 terms = [evalf(arg, prec + 10, options) for arg in v.args]
603 n = terms.count(S.ComplexInfinity)
604 if n >= 2:
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in <listcomp>(.0)
600 options['maxprec'] = min(oldmaxprec, 2*prec)
601
--> 602 terms = [evalf(arg, prec + 10, options) for arg in v.args]
603 n = terms.count(S.ComplexInfinity)
604 if n >= 2:
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_mul(v, prec, options)
648 from .numbers import Float
649 for arg in args:
--> 650 result = evalf(arg, prec, options)
651 if result is S.ComplexInfinity:
652 special.append(result)
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1494 try:
1495 rf = evalf_table[type(x)]
-> 1496 r = rf(x, prec, options)
1497 except KeyError:
1498 # Fall back to ordinary evalf if possible
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_exp(expr, prec, options)
892 def evalf_exp(expr: 'exp', prec: int, options: OPT_DICT) -> TMP_RES:
893 from .power import Pow
--> 894 return evalf_pow(Pow(S.Exp1, expr.exp, evaluate=False), prec, options)
895
896
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf_pow(v, prec, options)
833 # This determines the working precision that must be used
834 prec += 10
--> 835 result = evalf(exp, prec, options)
836 if result is S.ComplexInfinity:
837 return fnan, None, prec, None
~\Anaconda3\lib\site-packages\sympy\core\evalf.py in evalf(x, prec, options)
1491 Note that 0 is also represented as ``fzero = (0, 0, 0, 0)``.
1492 """
-> 1493 from sympy.functions.elementary.complexes import re as re_, im as im_
1494 try:
1495 rf = evalf_table[type(x)]
KeyboardInterrupt:
Evaluation of all terms to solve for the plot still takes too long. In Maple this takes quite a long time to solve as well.
Converting the equation to a numpy function speeds up the process
solution_in_frequency_domain_numeric = solution_in_frequency_domain.subs([(tau,2),(k,4),(m,1),(c,3/4)])
display(solution_in_frequency_domain_numeric)
\[\displaystyle \frac{i \left(1 - e^{4 i \pi \phi}\right) e^{- 8 i \pi \phi}}{2 \pi \phi \left(- 4 \pi^{2} \phi^{2} + 1.5 i \pi \phi + 4\right)}\]
solution_numeric = 2 * sp.integrate(sp.re(solution_in_frequency_domain_numeric*sp.exp(sp.I*2*sp.pi*phi*t)),(phi,0,10))
display(solution_numeric)
\[\displaystyle \frac{\int\limits_{0}^{10} \left(- \frac{\operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{- 4 \pi^{2} \phi^{2} e^{8 i \pi \phi} + 1.5 i \pi \phi e^{8 i \pi \phi} + 4 e^{8 i \pi \phi}}\right)}}{\phi}\right)\, d\phi + \int\limits_{0}^{10} \frac{\sin{\left(4 \pi \phi \right)} \operatorname{re}{\left(\frac{e^{2 i \pi \phi t}}{- 4 \pi^{2} \phi^{2} e^{8 i \pi \phi} + 1.5 i \pi \phi e^{8 i \pi \phi} + 4 e^{8 i \pi \phi}}\right)}}{\phi}\, d\phi + \int\limits_{0}^{10} \frac{\cos{\left(4 \pi \phi \right)} \operatorname{im}{\left(\frac{e^{2 i \pi \phi t}}{- 4 \pi^{2} \phi^{2} e^{8 i \pi \phi} + 1.5 i \pi \phi e^{8 i \pi \phi} + 4 e^{8 i \pi \phi}}\right)}}{\phi}\, d\phi}{\pi}\]
Analytical solution can be found using the residue theorem: https://groups.google.com/g/sympy/c/ZxgErmBe3qU. $\(\eqalign{ & \omega = {{\sqrt {247} } \over 8} + {3 \over 8}i \cr & \Delta \omega = {{\sqrt {247} } \over 4} \cr & \theta = \arctan \left( {{{{\mathop{\rm Im}\nolimits} \left( \omega \right)} \over {{\mathop{\rm Re}\nolimits} \left( \omega \right)}}} \right) \cr & P\left( \tau \right) = {{{e^{ - {\mathop{\rm Im}\nolimits} \left( \omega \right)\tau }}} \over {\left| \omega \right|}}\cos \left( {{\mathop{\rm Re}\nolimits} \left( \omega \right)\tau - \theta } \right) \cr & F\left( t \right) = {{ - 2} \over {\Delta \omega }}\left\{ {\matrix{ 0 & {t < 2} \cr {P\left( {t - 2} \right) - P\left( 0 \right)} & {2 \le t \le 4} \cr {P\left( {t - 2} \right) - P\left( {t - 4} \right)} & {t > 4} \cr } } \right. \cr} \)$
w_r = sp.sqrt(247)/8
w_i = 3/8
w_abs = sp.sqrt(w_r**2+w_i**2)
w_delta = sp.sqrt(247)/4
theta = sp.atan(w_i/w_r)
integral = sp.exp(-w_i*tau)/w_abs*(sp.cos(w_r*tau-theta))
solution_numeric_1 = -2/ w_delta * (integral.subs(tau,t-2)-integral.subs(tau,0))
solution_numeric_2 = -2/ w_delta * (integral.subs(tau,t-2)-integral.subs(tau,t-4))
p0 = sp.plotting.plot(solution_numeric_1, (t , 2 , 4),show=False)
p1 = sp.plotting.plot(solution_numeric_2, (t , 4 , 15),show=False)
p2 = sp.plotting.plot(0, (t , 0 , 2),show=False)
p0.append(p1[0])
p0.append(p2[0])
p0.show()
Manual integration (midpoint rule) is a lot faster:
import numpy as np
solution_numeric_eval = np.zeros(100)
length = 1.5 #upper limit phi
width = 0.01 #width of rectangles
t_eval = np.linspace(0,15,100)
solution_func = sp.lambdify((t,phi),2*sp.re(solution_in_frequency_domain_numeric*sp.exp(sp.I*2*sp.pi*phi*t)))
for i in np.linspace(width/2,length,int(length/width)):
solution_numeric_eval = np.add(solution_numeric_eval,solution_func(phi=i,t=t_eval)*width)
import matplotlib.pylab as plt
plt.plot(t_eval,solution_numeric_eval);
