Plotting the integral without getting negative domain results?

So I need to graph the magnitude of equation 2 (I labeled them in the comments) but I am not getting what is expected for the plots of the magnitudes of equation 2. It should be a graph that starts out at some y value that is high up on the y scale, and then over the frequency (which is the x value) gradually slope down in a curved fashion till it hits 0. an easy way to describe what I am supposed to get would be an upside down e^x graph. Does anyone have any suggestions as to what I need to do to my code to get these results? I have been struggling trying to get the right result now for a week. :/

if true
format long
a = 1 
b = 3E-7 
c = 5E-8 
f0 = 4E9 
sigma = 0.2
t0 = 0 
tmax = 2.*b 
f = linspace(1E6,1E10)
t = linspace(t0, tmax)
omega = 2.*pi.*f
omega0 = 2.*pi.*f0
yt = a.*exp((-(t-b).^2)./((2*c).^2))                       % equation (1)
fun = @(t) (yt.*exp(-1i.*omega.*t))
q = abs(integral(fun,0,6E-7, 'ArrayValued',1))             % equation (2)
figure(1)
plot(t,yt)                                                 % plot of (1) vs time for t:(0 to 6*10^-7)
figure(2)
plot(f,q)                                                  % plot of magnitude of F (eq 2) vs frequency for f:(10^6 to 10^10) 
r = q.^2                                                   % square of magnitude of F
figure(3)
plot(f,r)                                                  % plot of square of magnitude of F
q2 = integral(@(t) fun(t).^2,3.8E9,4.3E9, 'ArrayValued',1) % integral of square of magnitude of F
figure(4)
plot(f,q2)
  % code
end

Here's part of my code that includes what I was talking about above. Any input on what I could do to get the correct result will be highly appreciated!

Answer 1

[Walter Roberson]

A) You need to increase your resolution beyond the default number of points that linspace() generates

B) your fun is imaginary. squaring it is still imaginary. The integral of an imaginary quantity is imaginary and dependent upon the integration path. You should have noticed the

Warning: Imaginary parts of complex X and/or Y arguments ignored

[imarquez]

I put the first integral in absolute value form though, wouldn't that make the second integral also in absolute value form? I just updated my code to include it in the second integral anyways.

How do I increase the resolution?

[Walter Roberson]

Your second integral is

integral(@(t) fun(t).^2,3.8E9,4.3E9, 'ArrayValued',1)

which does not use the value of the first integral. And your fun(t) is imaginary.

Answer 2

[Star Strider]

See if substituting with this assignment does what you want:

q2 = integral(@(t) fun(t).*conj(fun(t)),3.8E9,4.3E9, 'ArrayValued',1);

[imarquez]

I had tried that earlier this morning and it just gave me a steeper form of the first graph that I plotted. Do you know how I can increase the resolution for the graphs? If I put it into logspace() my code doesn't work

[imarquez]

Actually never mind I figured out how to do the resolution. but I'm still puzzled as to why I'm not getting the graph that I'm expecting :/

[Star Strider]

What do you mean by ‘... an upside down e^x graph.’

Do you have an illustration?

[imarquez]

It would come out to something like the tachometer in this picture but think of the flat part starting on the left hand side of the graph or another way to think about it would be rotate the graph 180 degrees about the y axis but from the the center of the graph instead of on the y axis

http://www.sccs.swarthmore.edu/users/06/adem/engin/e58/lab3/index_files/image011.gif

[Star Strider]

That looks very much like a cumulative distribution function. It would be described as the integral of the integral I calculated in my Answer, with respect to the same independent variable.

I’m not certain how you would change your theory to accommodate the twice-iterated integral, but that’s how I would calculate it.

[imarquez]

I just don't want to do a second integration and have the output be incorrect for what we are analyzing. I did seem to get the graph I wanted now but I am running into another problem with an "Error using + Complex integer arithmetic is not supported" but only when I use it on both Ff and Zf and dont get the error when using it only on Zf.

This is what I have made my code into now and I get a graph that seems like it is correct when I put it in log scale, but I still cant get both graphs to run

if true
format long
a = 1 
b = 3E-7 
c = 5E-8 
f0 = 4E9 
sigma = 0.2
t0 = 0 
tmax = 2.*b 
f = linspace(1E6,1E10,1000)
t = linspace(t0, tmax,1000)
omega = 2.*pi.*f
omega0 = 2.*pi.*f0
G = -2.*pi.*f.*(pi.*c.^2.*f + sqrt(-1).*b)
D = (((b-(sqrt(-1)).*2.*pi.*c.^2.*f-t0)/(sqrt(2).*c)))
E = (((b-(sqrt(-1)).*2.*pi.*c.^2.*f-tmax)/(sqrt(2).*c)))
Ff = uint64(abs(sqrt(pi/2).*a.*c.*exp(G).*((-sqrt(-1).*(erfi(-sqrt(-1).*D)))-((-sqrt(-1).*(erfi(-sqrt(-1).*E)))))))       % equation (3)
Zf = uint64(((1-(sigma/2)).*Ff) + (sqrt(-1).*sqrt(pi/2).*((a.*c.*sigma)/4).*exp(-((2.*pi.*f+omega0).*(c.^2.*(2.*pi.*f+omega0)+2.*sqrt(-1).*b)))).*(-exp(4.*pi.*c.^2.*f.*omega0+2.*sqrt(-1).*b.*omega0).*((-sqrt(-1).*erfi(((tmax-b+sqrt(-1).*c.^2).*(2.*pi.*f-omega0))/(sqrt(2).*c)))-(-sqrt(-1).*erfi(((t0-b+sqrt(-1).*c.^2).*(2.*pi.*f-omega0))/(sqrt(2).*c))))+(-sqrt(-1).*(erfi(((tmax-b+sqrt(-1)*c.^2).*(2.*pi.*f+omega0))/(sqrt(-1).*c))))-(-sqrt(-1).*erfi(((t0-b+sqrt(-1).*c.^2).*(2.*pi.*f+omega0))/(sqrt(-1).*c)))))            % equation (5)
figure(6)                            
plot(f,Zf)                       % plot of Magnitude of G (eq 5) vs frequency for f:(10^6 to 10^10)
figure(7)
plot(f,Ff)
  % code
end

[Walter Roberson]

Addition is not supported between complex numbers and uint64().

It is possible to create uint64 values that are complex, but not by addition.

I have no idea why you are using uint64() in your code.

[imarquez]

In order to get a value out of Ff and Zf. Otherwise the numerical analysis will return NaN and 0 when they are just really small.

[Star Strider]

If double variables return NaN and 0, how do you know the uint64 values that do not are accurate? I would stay with double variables, and look closely at the code to see where the errors are if the results are not supposed to be 0. The minimum double-precision floating-point number is 2.2E-308.

[Walter Roberson]

So you are using uint64() just because uint64(0)/uint64(0) returns 0 instead of NaN? Why are you not using

X(isnan(X)) = 0;

replacing the NaN with 0

Or perhaps when you do the division A./B,

Bz = B;
Bz(Bz==0) = 1;
A./Bz

so that you divide by 1 instead of 0 in that location, getting 0 as the output ?