MATLAB Answers

Numerical integration in MATLAB

7 views (last 30 days)
Raj Patel
Raj Patel on 27 Oct 2020 at 17:46
Commented: Star Strider on 27 Oct 2020 at 21:33
Hello everyone,
I am trying to first numerically integrate a function and then plot the output of the integral. However, I have a variable T which is in the range (1:1000). I am not sure how to use a loop in integral and then plot values for all T with the corresponding value of integral. The MATLAB code is:
Hello everyone,
I am trying to first numerically integrate a function and then plot the output of the integral. However, I have a variable T in the integral which is in the range (1:1000). I am not sure how to use a loop in integral and then plot values for all T with the corresponding value of integral. 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;
T = 300;
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.

  1 Comment

Raj Patel
Raj Patel on 27 Oct 2020 at 18:07
The code is:
(Note: I removed a line where T was defined 2nd time)
Thanks in advance,
Raj Patel.
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)

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 27 Oct 2020 at 18:06
Try this:
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')
producing:
Since both ‘c’ and ‘fun’ are functions of ‘T’, they need to be coded as such.

  4 Comments

Show 1 older comment
Star Strider
Star Strider on 27 Oct 2020 at 18:44
As always, my pleasure!
It can be fully vectorised, avoiding the loop, using the 'ArrayValued' name-value pair:
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.
Raj Patel
Raj Patel on 27 Oct 2020 at 21:06
Oh! I see.
Thanks,
Raj Patel.
Star Strider
Star Strider on 27 Oct 2020 at 21:33
As always, my pleasure!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by