Need help plotting this function in Matlab

Question

I recently started using Matlab and I am trying to plot Imaginary part of a function. I can see I am making a mistake somewhere because I have the graph showing what I need to get and I am getting something else. Here is the link of the graph that I need to get, that I am getting and a function that I am plotting:

Wanted results, plotted function and obtained graph.

I know that this function has a singularity at frequency around 270 Hz and I don't know if a 'quadgk' can solve integral then. You guys probably know but theta function inside the integral is a heaviside and I don't know if that is causing problem maybe? Is there something wrong I am doing here? Here is the matlab code I wrote:

clear all
clc

ff=1:10:1000;

K1=(2*3)*log(2*cosh(135/6))/pi;
theta1=ones(size(ff));
theta1(ff<270)=0;

I1=zeros(size(ff));

for n = 1:numel(ff)
    f = ff(n);
    I1(n)= K1/f + (f/pi)*quadgk(@(x)(sinh(x/3)/(cosh(135/3)+cosh(x/3))-theta1(n))./((f^2)-4*(x.^2)), 0, inf);
end

plot(ff,I1, 'b');
hold on
axis([0 1000 -0.3 1])

I stared at the formula and your implementation in Matlab, there doesn't seem to be anything wrong with that (which the plot also suggests by the way). So I think you should try find an alternative for `quadgk` since that seems to be the problem. I can't help you with that though. — Fraukje, Sep 11 '13 at 13:27
I found this in the Matlab Documentation of `quad`(see http://www.mathworks.nl/help/matlab/ref/quad.html): "If the function is singular at points inside (a,b), write the integral as a sum of integrals over subintervals with the singular points as endpoints, compute them with `quadgk`, and add the results." However, it still seems as if there are multiple singularities in your function...? — Fraukje, Sep 11 '13 at 13:36
well, since `theta` is discontinuous, we'll find that step also in the `I1` function. — fpe, Sep 11 '13 at 13:43

Buck Thorn · Answer 1 · 2013-09-13T18:46:38.440

First, I altered the expression for your heaviside function to what I think is the correct form.
Second, I introduced the definitions of mu and T explicitly.
Third, I broke down the integral into 4 components, as follows, and evaluated them individually (along the lines of what Fraukje proposed)
Fourth, I use quadl, although this is prob. of minor importance.
(edit) Changed the range of ff

Here's the code with changes:

fstep=1;
ff=[1:fstep:265 275:fstep:1000];

T = 3;
mu = 135;

df = 0.01;
xmax = 10;

K1=(2*T/pi)*log(2*cosh(mu/(2*T)));
theta1=ones(size(ff));
theta1(ff-2*mu<0)=0;

I1=zeros(size(ff));    
for n = 1:numel(ff)
    f = ff(n);
    sigm1 = @(x)  sinh(x/T)./((f^2-4*x.^2).*(cosh(mu/T)+cosh(x/T)));
    sigm2 = @(x)   -theta1(n)./(f^2-4*x.^2);
    I1(n) = K1/f  + (f/pi)*quadl(sigm1,0,f/2-df);      % term #1
%    I1(n) = I1(n) + (f/pi)*quadl(sigm1,f/2+df,xmax);   % term #2
%    I1(n) = I1(n) + (f/pi)*quadl(sigm2,0,f/2-df);     % term #3
%    I1(n) = I1(n) + (f/pi)*quadl(sigm2,f/2+df,xmax);   % term #4
end

I selected to split the integrals around x=f/2 since there is clearly a singularity there (division by 0). But additional problems occur for terms #2 and #4, that is when the integrals are evaluated for x>f/2, presumably due to all of the trigonometric terms.

If you keep only terms 1 and 3 you get something reasonably similar to the plot you show:

enter image description here

However you should probably inspect your function more closely and see what can be done to evaluate the integrals for x>f/2.

EDIT

I inspected the code again and redefined the auxiliary integrals:

I1=zeros(size(ff));
I2=zeros(size(ff));
I3=zeros(size(ff));

for n = 1:numel(ff)
    f = ff(n);
    sigm3 = @(x)  sinh(x/T)./((f^2-4*x.^2).*(cosh(mu/T)+cosh(x/T))) -theta1(n)./(f^2-4*x.^2);
    I1(n) = K1/f  + (f/pi)*quadl(sigm3,0,f/2-df); 
    I2(n) = (f/pi)*quadl(sigm3,f/2+df,10);
end
I3=I2;
I3(isnan(I3)) = 0;
I3 = I3 + I1;

This is how the output looks like now:

enter image description here

The green line is the integral of the function over 0<x<f/2 and seems ok. The red line is the integral over Inf>x>f/2 and clearly fails around f=270. The blue curve is the sum (the total integral) excluding the NaN contribution when integrating over Inf>x>f/2.

My conclusion is that there might be something wrong with the curve as you expect it to look.

score 1 · Answer 2 · answered Sep 11 '13 at 14:03

So far I'd proceed this way:

clc,clear

T = 3;
mu = 135;

f = 1E-04:.1:1000;

theta = ones(size(f));
theta(f < 270)= 0;

integrative = zeros(size(f));
for ii = 1:numel(f)
    ff = @(x) int_y(x, f(ii), theta(ii));
    integrative(ii) = quad(ff,0,2000);
end

Imm = ((2*T)./(pi*f)).*log(2*cosh(mu/(2*T))) + (f/pi).*integrative;

Imm1 = exp(interp1(log(f(1:2399)),log(Imm(1:2399)),log(f(2400):.001:f(2700)),'linear','extrap'));
Imm2 = exp(interp1(log(f(2985:end)),log(Imm(2985:end)),log(f(2701):.001:f(2984)),'linear','extrap'));


plot([(f(2400):.001:f(2700)) (f(2701):.001:f(2984))],[Imm1 Imm2])
hold on
axis([0 1000 -1.0 1])
plot(f,Imm,'g')
grid on
hold off

with

function rrr = int_y(x,f,theta)
T = 3;
mu = 135;
rrr = ( (sinh(x./T)./(cosh(mu/T) + cosh(x/T))) - theta ) ./ (f.^2 - 4.*(x.^2));
end

I've come up with this plot:

enter image description here

Need help plotting this function in Matlab

2 Answers2