|
Size: 3313
Comment: Add Cauchy Residue Theorem
|
← Revision 3 as of 2021-06-24 22:30:22 ⇥
Size: 3342
Comment: forgot factor 2 pi i
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 75: | Line 75: |
| pretty_print(HtmlFragment(r' |
pretty_print(HtmlFragment(r'$\int_\gamma f(z)\: d z = 2\pi i \left(%s\right)$'% |
| Line 84: | Line 84: |
| {{attachment:geodesics1.png}} | {{attachment:residues_lemniscate.png}} |
Cauchy residue theorem
Pablo Angulo. Shows a function f of a complex variable, a plane curve gamma given as a map R -> C t -> gamma(t) and plots the contribution of the residues of f to the contour integral of f along gamma
xxxxxxxxxximport scipy.integrate as integfrom sage.misc.html import HtmlFragmentvar('z t')def winding_number(gamma, z0): t, = gamma.variables() assume(t, 'real') x,y = gamma.real_part(), gamma.imag_part() z0 += I*0 x0,y0 = z0.real_part(), z0.imag_part() xmx0 = x - x0 ymy0 = y - y0 r2 = xmx0^2 + ymy0^2 xp = xmx0.derivative(t) yp = ymy0.derivative(t) integrando(t) = (xmx0*yp -ymy0*xp)/r2 i,e = numerical_integral(integrando, 0, 2*pi) return round(i/(2*pi))def cauchy_residue(f, gamma): z, = f.variables() t, = gamma.variables() ftheta = f.subs({z:gamma})*gamma.diff(t) quad_integral = ( integ.quad(lambda t0:ftheta(t=t0).real_part(), 0, 2*pi)[0], integ.quad(lambda t0:ftheta(t=t0).imag_part(), 0, 2*pi)[0], ) poles = [ d[z] for d in solve(1/f==0,z, solution_dict=True) ] residues = [ #pole, residue and its winding number wrt gamma (z0, f.residue(z==z0), winding_number(gamma, z0)) for z0 in poles ] residue_integral = 2*pi*i*sum( res*win for z0,res,win in residues ).simplify_full() return residues, residue_integral, quad_integraldef cauchy_residue_plot(f, gamma): residues, residue_integral, quad_integral = cauchy_residue(f, gamma) max_win = max(win for z0,res,win in residues) min_win = min(win for z0,res,win in residues) return ( complex_plot(f,(-3,3),(-3,3)) + point2d([(z0.real_part(), z0.imag_part()) for z0,res,win in residues], size=50, zorder=5, color='white') + sum([point2d((z0.real_part(), z0.imag_part()), size=30, zorder=6, color=(max(0,win/(max_win-min_win)),max(0,-win/(max_win-min_win)),0)) for z0,res,win in residues]) + parametric_plot( (gamma.real_part(), gamma.imag_part()), (t, 0, 2*pi)) + sum([text('$%s$'%latex(res),(z0.real_part()+0.1, z0.imag_part()+0.3), fontsize=10, color='white', background_color=(max(0,win/(max_win-min_win)),max(0,-win/(max_win-min_win)),0), alpha=1, zorder=4) for z0,res,win in residues]) )def _(f = cos(z)/(z^4-1), #Lemniscate gamma = 2*(cos(t)/(1+sin(t)^2) + i*sin(t)*cos(t)/(1+sin(t)^2))): residues, residue_integral, quad_integral = cauchy_residue(f, gamma) show('Residues:', residues) pretty_print(HtmlFragment(r'$\int_\gamma f(z)\: d z = 2\pi i \left(%s\right)$'% (' + '.join(r'%s\cdot %s'%(latex(win), latex(res)) for z0,res,win in residues)))) pretty_print(HtmlFragment(r'$\int_\gamma f(z)\: d z = %s$'% residue_integral.simplify_full())) pretty_print(HtmlFragment(r'$\int_\gamma f(z)\: d z \approx %.8f + %.8f*I$'% quad_integral)) cauchy_residue_plot(f, gamma).show(figsize=(8,8))