Plot integral of a function
显示 更早的评论
Hello,
I'm tryign to plot the integral of the function rho in the following code, but I'm having trouble with plotting it. Perhaps anyone can help me with the code for plotting the integral of this rho function in the different areas? This is what I've tried:
q = 1.5*10^-16;
a = 1;
Na = 10^16;
xn = 3;
xp = -2;
rho_righ = @(x) a.*x;
rho_left = @(x) -q*Na;
nop = @(x) 0;
figure(1);
fplot(rho_righ, [0,xn]);
hold on;
fplot(rho_left, [xp,0]);
hold on;
fplot(nop, [-5,xp]);
hold on;
fplot(nop, [xn,5]);
grid on;
legend('{\rho}=qax','{\rho}=-qNa');
hold off;
upper_limit = linspace(-xp, 0);
upper_limit2 = linspace(0, xn);
xval = arrayfun(@(uplim) integral(rho, 0, uplim, 'ArrayValued',true), upper_limit);
xval2 = arrayfun(@(uplim) integral(rho, 0, uplim, 'ArrayValued',true), upper_limit2);
figure(2)
plot(upper_limit, xval);
hold on;
plot(upper_limit2, xval2);
grid on;
Thank you very much!
采纳的回答
Star Strider
2020-12-27
The problem appears to be that ‘rho’ itself does not exist.
I have no idea what you want ‘rho’ to be, however defining it as:
rho = @(x) rho_left(x) .* (x <= 0) + rho_righ(x) .* (x > 0);
will likely be an appropriate place to begin, and that produces figure(2) (and myriad Warning messages).. Whatever you intend it to be, it needs to be defined.
4 个评论
Elinor Ginzburg
2020-12-28
Star Strider
2020-12-28
编辑：Star Strider
2020-12-28
My pleasure!
First, arrayfun has its uses, however it’s simply a for loop in disguise, so it’s more efficient to just use a for loop. Also, the ‘&&’ is a ‘short-circuit’ operator for strictly logical values. For vectors, use a single ‘&’ operator.
Second, I’m not certain what you’re doing in the second integration. It looks as though you’re doing a double integration over the same limits, where ‘V’ is integrating ‘EF’ is integrating ‘rho’. There might be ways to make that more tractable if I understand what you’re doing.
Meanwhile, while I wait for you to clarify what the the ‘voltage’ integration is (possibly to make it faster and more efficient), your (slightly-revised) code is:
q = 3;
a = 1/2;
Na = 1;
xn = 2;
xp = -1;
rho_righ = @(x) q*a.*x;
rho_left = @(x) -q*Na;
rho = @(x) rho_left(x) .* ((x <= 0) & (x > xp)) + rho_righ(x) .* ((x > 0) & (x <= xn));
figure(1);
fplot(rho);
grid on;
xlabel('x [{cm^{-3}}]');
ylabel('{\rho} [{cb/cm^{-3}}]');
title('Charge Density Distribution');
upper_limit = linspace(-10, 10);
EF = @(uplim) integral(rho, 0, uplim, 'ArrayValued',true) - 3;
for k = 1:numel(upper_limit)
Electrical_field(k) = EF(upper_limit(k));
end
figure(2)
plot(upper_limit, Electrical_field);
grid on;
xlabel('x [{cm^{-3}}]');
ylabel('E [V/cm]');
title('Electrical Field');
V = @(uplim) integral(EF, 0, uplim, 'ArrayValued',true);
for k = 1:numel(upper_limit)
voltage(k) = V(upper_limit(k));
end
figure(3)
plot(upper_limit, voltage);
grid on;
xlabel('x [{cm^{-3}}]');
ylabel('V [V]');
title('Voltage');
EDIT — (28 Dec 2020 at 13:28)
It only ended up taking 1647.612439 seconds (27.46 minutes) for all 100 loop iterations. In the end, it produced:
To save you the trouble of integrating it, the ‘voltage’ vector is:
voltage = [1.49999729339940 1.49999696906513 1.50015522691315 1.50000236066725 1.49999973680205 1.50000048145849 1.49999931842135 1.49999972895611 1.49999906076886 1.49999918636686 1.49999895100276 1.49999918002933 1.49999971120408 1.49999952794135 1.49999962815786 1.50000030206563 1.49999982037269 1.50000168699648 1.50000016177641 1.49999957361846 1.49999976226557 1.49999946542346 1.50000107234303 1.49999866049175 1.50000007607924 1.49999936777078 1.49999950297102 1.50000078110427 1.49999943944442 1.50000006989847 1.50000061270948 1.49999949173887 1.49999992644981 1.49999988813672 1.49999975368376 1.49999941182790 1.49999970577751 1.50000030183291 1.49999923594608 1.49999958590485 1.49999927394898 1.49999985898082 1.50000059556414 1.49999973974386 1.50000074795001 1.48760330578512 1.37128864401592 1.13253749617386 0.771349862258953 0.287725742271197 -0.302772650492271 -0.902134290564043 -1.48294494789750 -2.03283730066711 -2.53944402704733 -2.99039780521262 -3.37333131333744 -3.67587722959624 -3.88566823216349 -3.99033699921365 -4.00000064537387 -4.00000095524024 -3.99999924718646 -4.00001441974673 -3.99999995888804 -4.00000125703370 -4.00000118209333 -4.00000010402610 -3.99999897834425 -4.00000301314983 -3.99999985682291 -3.99999997325554 -4.00000049834742 -3.99999985595181 -4.00000027641440 -4.00000226576712 -4.00000150165796 -4.00000056720757 -3.99999992990422 -4.00000016104336 -3.99999930707085 -3.99999881231945 -4.00000009175432 -3.99999780633590 -3.99999761823918 -4.00000016397481 -3.99999902244021 -3.99999924718711 -3.99999907700015 -4.00000099591624 -3.99999769551857 -4.00000009842234 -3.99999631419898 -3.99999899595210 -3.99999286143871 -3.99999971820702 -3.99999724890175 -3.99999932916977 -3.99999697121052 -3.99999797250258];
I was able to recover it from my Workspace window and copy it to here.
Elinor Ginzburg
2020-12-28
Star Strider
2020-12-28
As always, my pleasure!
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 MATLAB 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
