Numerical Integration In Matlab

I need some help solving numerical integration.

The integral is:

fun = @(x) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* 1 .* x)).^2 ./ sqrt(1 - (0.0003 .* 1 .* x).^2))); Note: asind = sin^-1(....) --> (I am not sure if this is correct)

C = @(T) integral(fun, 0, 517.39 ./ T);

The limit of x is from: x = 0 to x = (500/T)

T ranges from (1 - 1000) (which is T = logspace(0, 3))

Can anyone help me solve this integral. I want to plot a loglog plot of (C as a function of T). I am getting some errors which I am not able to overcome.

Thanks in advance.

Raj.

Answer 1

[Star Strider]

The integral function will produce one value for the specific limits of integration. To create a vector from its outputs, a loop is necessary.

Try this:

fun = @(x) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* 1 .* x)).^2 ./ sqrt(1 - (0.0003 .* 1 .* x).^2)));
T = logspace(0, 3);
C = @(T) integral(fun, 0, 517.39 ./ T);

for k = 1:numel(T)
    Cint(k) = C(T(k));
end

figure
plot(T, Cint)
grid

.

[Raj Patel]

Oh I am sorry, the function is:

fun = @(x) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* T .* x)).^2 ./ sqrt(1 - (0.0003 .* T .* x).^2)));

Can you see the changes?

Thanks

Raj

[Star Strider]

My pleasure!

No worries!

The ‘fun’ function now needs to be coded as:

fun = @(x,T) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* T .* x)).^2 ./ sqrt(1 - (0.0003 .* T .* x).^2)));

and ‘C’ now needs to be made compatible with it:

C = @(T) integral(@(x)fun(x,T), 0, 517.39 ./ T);

Then, with those changes, the complete code is now:

fun = @(x,T) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* T .* x)).^2 ./ sqrt(1 - (0.0003 .* T .* x).^2)));
T = logspace(0, 3);
C = @(T) integral(@(x)fun(x,T), 0, 517.39 ./ T);

for k = 1:numel(T)
    Cint(k) = C(T(k));
end

figure
plot(T, Cint)                                   % Linear ‘x’
grid

figure
semilogx(T, Cint)                               % Logarithmic ‘x’
grid

I forgot about the possibility of a logarithmic scale for the x-axis previously. I include it here as an option. .

[Raj Patel]

Thank you so much. Apprecite your effort.

Raj.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Raj Patel]

I wanted to ask that if I have a prefactor of (5 * (1/T)) before the integral, where should I put it? I tried putting it like this, but it is not working.

fun = @(x,T) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* T .* x)).^2 ./ sqrt(1 - (0.0003 .* T .* x).^2)));

T = logspace(0, 3);

C = 2.35 .* 10.^-23 .* (1 ./ T) .* @(T) integral(@(x)fun(x,T), 0, 517.39 ./ T);

Thanks,

Raj.

[Star Strider]

As always, my pleasure!

The code to do that needs to be something like this:

fun = @(x,T) (x .* (exp(x) ./ (exp(x) - 1).^2) .* ((asind(0.0003 .* T .* x)).^2 ./ sqrt(1 - (0.0003 .* T .* x).^2)));
T = logspace(0, 3);
C = @(T) 2.35 .* 1E-23 .* (1 ./ T) .* integral(@(x)fun(x,T), 0, 517.39 ./ T);
for k = 1:numel(T)
    Cint(k) = C(T(k));
end

figure
semilogx(T, Cint)                                  % Logarithmic ‘x’
grid

I changed the ‘C’ function to be compatible with the premultiplication requirement. Note that now ‘C’ is a second anonymous function. An explicit for loop is still the most efficient option here.

[Raj Patel]

Thanks Star.

Raj

[Star Strider]

As always, my pleasure!

[Raj Patel]

Hi Star,

I wanted some help with a Matlab coding. I am trying to integrate a function and this function has a variable T which ranges from 1:1000. How do I get numerical integral for all these values of T and then plot it as a function of T? Could you help me? The Matlab code is:

kb = 1.38 .* 10.^-23;
h = 1.05 .* 10.^-34;
wd = 6.94 .* 10.^13;
vs = 6084;
B1 = 2.7 .* 10.^-19;
B2 = 170;
A1 = 2 .* 10.^-45;
D = 0.004;
c = (h.^2 .* D) ./ (2 .* pi * pi * vs * kb .* T .* T);
T = 1:1000;
fun = @(x) ((x.^4 .* exp((h.*x) ./ (kb .* T))) ./ (((exp((h.*x) ./ (kb .* T))-1).^2) .* ((D.*B1 .* x.^2 .* T .* exp(-B2/T)) + (D.*A1 .* x.^4) + vs)));
K = c .* integral(fun,0,wd)
plot(T, K)

Thanks in advance,

Raj Patel.

[Star Strider]

As always, my pleasure!

I thought we just solvedd that exact problem? I only posted the few lines that changed.

The full code is:

kb = 1.38 .* 10.^-23;
h = 1.05 .* 10.^-34;
wd = 6.94 .* 10.^13;
vs = 6084;
B1 = 2.7 .* 10.^-19;
B2 = 170;
A1 = 2 .* 10.^-45;
D = 0.004;
T = 300;
c = @(T) (h.^2 .* D) ./ (2 .* pi * pi * vs * kb .* T .* T);
fun = @(x,T) ((x.^4 .* exp((h.*x) ./ (kb .* T))) ./ (((exp((h.*x) ./ (kb .* T))-1).^2) .* ((D.*B1 .* x.^2 .* T .* exp(-B2/T)) + (D.*A1 .* x.^4) + vs)));

T = 1:1000;
for k = 1:numel(T)
    K(k) = c(T(k)) .* integral(@(x)fun(x,T(k)),0,wd);
end

figure
plot(T, K)
grid
xlabel('T')
ylabel('K')

It produces the same plot, so I will not re-post it.

This also works (without the loop, using 'ArrayValued'):

c = @(T) (h.^2 .* D) ./ (2 .* pi * pi * vs * kb .* T .* T);
fun = @(x,T) ((x.^4 .* exp((h.*x) ./ (kb .* T))) ./ (((exp((h.*x) ./ (kb .* T))-1).^2) .* ((D.*B1 .* x.^2 .* T .* exp(-B2./T)) + (D.*A1 .* x.^4) + vs)));

T = 1:1000;

K = c(T) .* integral(@(x)fun(x,T),0,wd,'ArrayValued',1);

figure
plot(T, K)
grid
xlabel('T')
ylabel('K')

producing the same plot.

The ‘fun’ code changed slightly to make it fully vectorised.

.