How can I code a contour integral representation of cosine?
8 views (last 30 days)
% "Contour Integral Representation of Cosine"
% making 'z' a possible number such as pi
N = 10000000; % example number
z_lower = 0;
z_upper = 6*pi;
z1 = linspace(z_lower,z_upper,N);
y = 1;
Sl = y - 10*1i;
Su = y + 10*1i;
ds = (Su - Sl)/N; % size of each step
sum = 0.0;
for z = linspace(z_lower,z_upper,N)
for Sn = linspace(Sl,Su,N)
sum = sum + ((exp(Sn) - (z.^2/4*Sn))/sqrt(Sn))*ds;
sum = sum*(sqrt(pi)/2*pi*1i);
The above is my poor attempt! I am trying to write a code to numerically approximate the integral along a suitable contour and plot the approximation against the exact value of cos z. I would like to do this twice - once for z ∈ [0, 6π] and once for complex valued z in the range z ∈ [0 + i, 6π + i]. Then I intend to plot both the real and imaginary part of the computed cos z.
I have a following lost of steps that I am trying to implement!
- Choose γ, Sl, Su , N.
- Step through z from z lower to z upper (using a different number of steps (other than N) for this).
- For each value of z compute cos z by stepping along the contour in N discrete steps from Sl to Su .
- For each value of sn along the contour compute the integrand "((exp(Sn) - (z.^2/4*Sn))/sqrt(Sn))" [I have attached an image to make this integrand clearer to read] and add it to the rolling sum.
- Plot the result.
- Repeat using a better integration rule
Image of the integrand!
I am pretty new to general forms of computation and would really appreciate any help! Thanks!
David Goodmanson on 24 Oct 2022
Edited: David Goodmanson on 24 Oct 2022
this integral converges very slowly along the vertical path gamma-i*inf to gamma+i*inf.
Convergence depends on what happens for large values of s. In that case, you can drop the z^2/s term in the exponential and look at exp(s). When the path of integration has large imaginary part, exp(s) is oscillatory and only 1/sqrt(s) is making the integrand smaller. Convergence is slow.
Rather than stick to the vertical path, which will prove the point philosophically but is difficult to achieve numerically, it helps to use a different one. The path is moved over from the original path in a continuous fashion. Here the new path is a rectangle:
o ^ o = origin
The idea is that at the left hand parts of the path, s has a large negative real part and exp(s) dies off quickly rather than oscillating.
In the code below it would have been nice to use the set
I1 = integral(@(t) fun(t,z), -inf-i*a, b-i*a);
I2 = integral(@(t) fun(t,z), b-i*a, b+i*a);
I3 = integral(@(t) fun(t,z), b+i*a, -inf+i*a);
but the integral function can't handle complex end points at infinity, so the vertical offset is inserted manually in I1 and I3. With the choices for a and b the code does well up to z = 7*pi. You can mess around with a and b to see what works.
z = pi/3;
a = 3
b = 1
I1 = integral(@(t) fun(t,z,-i*a), -inf, b);
I2 = integral(@(t) fun(t,z,0 ), b-i*a, b+i*a);
I3 = integral(@(t) fun(t,z, i*a), b, -inf);
I = sqrt(pi)/(2*pi*i)*(I1+I2+I3)
delta = I - cos(z) % check
function y = fun(t,z,u)
s = t+u;
y = (exp(s-z.^2./(4*s))./sqrt(s));
I = 0.5000 + 0.0000i
delta = -2.7756e-16 + 6.2638e-17i