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!

Answer 1

[Star Strider]

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.

[Elinor Ginzburg]

Thank you very much for your answer, I totally missed that. I managed to plot the first integral of rho, and now I'm also trying to plot the second integral of this, but for some reason it starts to run and doesn't finish even after few minutes. perhaps you have any idea what could I do wrong? here's my new code:

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;
Electrical_field = arrayfun(EF, upper_limit);

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);
voltage = arrayfun(V, upper_limit);
figure(3)
plot(upper_limit, voltage);
grid on;
xlabel('x [{cm^{-3}}]');
ylabel('V [V]');
title('Voltage');

Thank you very much for your help!

[Star Strider]

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]

Thank you so much!

yeah that was exactly what I ment, and the results match my calculations.

Thank you very much for your help and advices!

[Star Strider]

As always, my pleasure!