How can I code a contour integral representation of cosine?

% "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!

[Torsten]

Could you explain what this has to do with a contour integral representation of cos(z) ? I don't get it.

And what is z in your code ? z1 ?

And as written, you only get a single value for sum. I think this is not what you want.

Maybe you could explain better what ((exp(Sn) - (z.^2/4*Sn))/sqrt(Sn)) has to do with cos(z).

[Richard Mc Enteggart]

Thanks for getting back to me! I hope I can explain things a bit clearer.

I have inserted an image of a cosine function plot. I would like to numerically approximate the integral along a suitable contour and plot the approximation against the exact value of cos(z). There is nothing special with the letter "z" it can be any letter in particular. By using a McLaurin expansion on paper we could get something like cos z = 1 − (z^2)/2 + (z^4)/4! +.....

However I am looking to write cos(z) in terms of an integral in the complex plane which I have supplied here with fig. 1.

Essentially, I am choosing an equi-sampled set of points sn along the contour, from Sl to Su and separated by ∆s. This is Euler’s method.

The image of the graph I have supplied is what I would be expecting-ish.

Figure 1 - 'cos(z) written as an integral in the complex plane'

[Torsten]

Here is a test for your formula to experiment with. It seems necessary to choose large values for N to get a good approximation.

gamma = 0.1;
z = pi;
N = 50;
f = @(t)-1i/(2*sqrt(pi))*exp((gamma+1i*t)-z^2./(4*(gamma+1i*t))) * 1i ./ sqrt(gamma+1i*t);
sol = integral(f,-N,N)

sol = -1.0747 - 0.0000i

[Richard Mc Enteggart]

Hiya! thank you for getting back to me. I have ran that lil segment of code you just send to me and when N is much larger, it does seem to change (probably as a better approximation!) while imaginary is always zero it seems.

So I have been headscratching this problem for a while now! I have completed the contour section. This code just below will give me the exact values for cos(z) over one full period, so this is my contour or "reference" in order to compare to the approximation!

clear, clc, close all

% Work with one full period

z = 0 : .1 : 2*pi;

y = cos(z);

% Plot the goal

plot(z, y, 'g', 'Linewidth', 3)

title('Study of Maclaurin series')

xlabel('z')

ylabel('cos(z) with different number of terms')

axis([-2*pi 2*pi -3 3])

grid on

hold on

% Consider 4 terms

smp = maclaurin_cos(z, 4);

plot(x, smp, 'ro')

% Consider 6 terms

smp = maclaurin_sin(z, 6);

plot(z, smp, 'b-.')

% Consider 10 terms

smp = maclaurin_cos(z, 10);

plot(z, smp, 'k', 'linewidth', 3)

% Label the calculated lines

legend('cos(z)', '4-term series', ...

'6 terms', '10 terms')

This is the waveform snippet I have of the exact values of cos(z) over one full period

If you remember, I am trying to get this!

I already have "exact" here, I am trying to get the "real" and "imaginary" waveforms (which I think that snippet of code you submitted could be used!)

The approximation I want to compute is here - this is the integrand and how it relates to cos(z). I want to compute this integrand (and plot it!), which is the approximation and plot it in MATLAB. This part is what is giving me the massive headaches!

[Torsten]

If you need to approximate the integral with midpoint rule, I leave this to you.

This is what you should get - I evaluated the integral directly with gamma = 0.1 in the interval [gamma-50*i:gamma+50*i].

gamma = 0.1;
z = 0:0.1:2*pi;
N = 50;
f = @(t,z)-1i/(2*sqrt(pi))*exp((gamma+1i*t)-z^2./(4*(gamma+1i*t))) * 1i ./ sqrt(gamma+1i*t);
cos_approx = arrayfun(@(z)integral(@(t)f(t,z),-N,N),z);
hold on
plot(z,real(cos_approx))
plot(z,imag(cos_approx))
plot(z,cos(z))
hold off

[Richard Mc Enteggart]

Wow! That is super interesting. Could you tell me how you implemented this? I copied and pasted it into MATLAB, it runs but it did not give me any waveforms! I simply opened a new script, added that block of code you just gave me and it seems to run as there are no errors in the console or variable sections, but it does not give me a plot. (No warnings either! Its showing me a green tick on the top right).

[Torsten]

I don't know. As you can see, it runs under MATLAB R2022b. Which MATLAB version do you have in use ?

[Richard Mc Enteggart]

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

Answer 1

[David Goodmanson]

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

[Richard Mc Enteggart]

Thanks for getting back to me! That is incredibly interesting! Your expalantions are making sense. I did not know my integrand was getting significantly smaller.

Is there any way I could plot this bloc of code?

[David Goodmanson]

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.