# Can someone help me check if I did this correctly?

1 view (last 30 days)
on 1 Dec 2023
Answered: Torsten on 1 Dec 2023
Here's the problem:
And the script:
syms r a0;
% Part (a) - Probability Density Function p(r) and its integral from 0 to infinity
p_r = 4 / a0^3 * r^2 * exp(-2 * r / a0);
integral_p = int(p_r, r, 0, Inf);
disp("Part (a) - Integral of p(r) from 0 to infinity:");
Part (a) - Integral of p(r) from 0 to infinity:
disp(integral_p);
% Part (b) - Limit as r approaches infinity
limit_r_infinity = limit(p_r * r^2 / exp(2 * r / a0), r, Inf);
disp("Part (b) - Limit as r approaches infinity:");
Part (b) - Limit as r approaches infinity:
disp(limit_r_infinity);
% Part (c) - Visualization of p(r) and its maximum
a0 = 5.59e-11;
p_r = @(r) 4 / a0^3 * r.^2 .* exp(-2 * r / a0);
r_vals = linspace(0, 2*a0, 1000);
p_values = p_r(r_vals);
[max_value, max_index] = max(p_values);
max_r = r_vals(max_index);
plot(r_vals, p_values);
hold on;
plot(max_r, max_value, 'ro');
text(max_r, max_value, sprintf('(%.2e, %.2e)', max_r, max_value), 'VerticalAlignment', 'bottom');
title('Probability Density Function p(r)');
xlabel('r');
ylabel('p(r)');
legend('p(r)', 'Maximum');
hold off;
% Part (d) - Calculation of P(4a0)
P_4a0 = integral(p_r, 0, 4 * a0);
disp("Part (d) - Integral P(4a0):");
Part (d) - Integral P(4a0):
disp(P_4a0);
0.9862
% Part (e) - Expected value mu
syms r a0;
% Define the integrand for the expected value mu
integrand_mu = r^3 * exp(-2 * r / a0);
% Calculate the expected value mu
mu = 4 / a0^3 * int(integrand_mu, r, 0, Inf);
disp("Part (e) - Expected value mu:");
Part (e) - Expected value mu:
disp(mu);

Torsten on 1 Dec 2023
Replace
syms r a0
by
syms r a0 positive
to get the limits resolved.