How can I code a contour integral representation of cosine?

9 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;
end
end
sum = sum*(sqrt(pi)/2*pi*1i);
plot(Sn,sum)
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!
  8 Comments
Richard Mc Enteggart
Richard Mc Enteggart on 24 Oct 2022
Edited: Richard Mc Enteggart on 24 Oct 2022
I use MATLAB R2022b as well! I ran the code suggested by you and it works! My machine just did not compute it I guess haha
Thank you very much :D

Sign in to comment.

Answers (1)

David Goodmanson
David Goodmanson on 24 Oct 2022
Edited: David Goodmanson on 24 Oct 2022
Hi Richard,
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:
-inf+ia a+ib
------<-------
o ^ o = origin
------>------^
-inf-ia a-ib
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));
end
I = 0.5000 + 0.0000i
delta = -2.7756e-16 + 6.2638e-17i
  2 Comments
David Goodmanson
David Goodmanson on 24 Oct 2022
Edited: David Goodmanson on 25 Oct 2022
Hi Richard, I expanded the original answer to hopefully provide some clarification. As to plotting, if you use the 'arrayvalued' option in the integral function then z can be a vector of values and you can plot the resulting output vs z.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by