problem with integrating a piece-wise function

Hello everybody,

I've been trying to integrate a piece-wise function using cumsum and for loop, but this doesn't work correctly. Mathematically/Physically, the problem is defined as follows: I have a charge density rho(x) vs. coordinate x between x_min and x_max across 2 materials with different dielectric permittivities eps1 and eps2.

For a uniform media (1 material) I find an electric field E by summing up rho(x)/eps over small intervals dx:

x =x_min:dx:x_max;
rho = @(x) rho_A*[heaviside(x-x_min)-heaviside(x-x_max)]
E=cumsum(rho(x(k))/eps_r)*dx;

That works fine. However, with 2 materials when I use for loop and if else statement to assign different eps to the 2 media, namely:

x =x_min:dx:x_max;
for k = 1:length(x)    
    if x(k)<=0;
        eps_r=10;
        rho = @(x) rho_A*(heaviside(x-x_min)-heaviside(x-x_max));
        E(k)=cumsum(rho(x(k))/eps_r)*dx;
    else
        eps_r=15;
        rho = @(x) rho_A*(heaviside(x-x_min)-heaviside(x-x_max));
        E(k)=cumsum(rho(x(k))/eps_r)*dx;
    end
end

That doesn't work correctly - cumsum doesn't sums up rho(x)/eps*dx over x =x_min:dx:x_max, but instead mimics the charge density distribution rho(x). I guess I am not using the right tools for the job. Please help to solve the problem.

Many thanks

Answer 1

[Star Strider]

You’re not assigning different values for ‘eps’:

        eps_r=15;                                                   % <- ‘eps_r’ Here
        rho = @(x) rho_A*(heaviside(x-x_min)-heaviside(x-x_max));
        E(k)=cumsum(rho(x(k))/eps)*dx;                              % <- ‘eps’ Here

Also, please do not use ‘eps’ as a variable name. It is the name of a built-in MATLAB function, so if you want to use it as a function later in your code, you will not be able to. This is called ‘overshadowing’, and it is best not to do it. The solution is to rename your ‘eps’ variable to ‘epsln’ or something that makes sense in the context of your problem.

[Anton Velychko]

sorry, I made a mistake in the code of the original post when editing it. I've corrected it now. The actual code I use doesn't have this mistake and using different epsilon (eps and eps_r) is not the cause of the problem. The problem of not summing up rho(x)/eps*dx is still there.

[Star Strider]

I think I see the problem. You’re not actually summing over the ‘x’ vector, but computing the value of ‘E’ at each value of ‘x’. The cumsum call is summing one element, so it is having no effect, since I believe all your other arguments to it are also scalars.

I’m not certain what you’re doing, and I don’t know what your other variables should be, so I made some up. This code works as I believe you want it to, but I have no idea if it produces the correct answers. (I would use cumtrapz rather than cumsum, so I put a commented line in for it in the event you want to experiment with it.)

The code:

x_min = -10;
x_max =  10;
dx    = 0.1;
rho_A =  pi;
rho = @(x) rho_A*(heaviside(x-x_min)-heaviside(x-x_max));
x = x_min:dx:x_max;
% E = @(x,eps_val) [(cumtrapz(rho(x)/eps_val(1))*dx).*(x <= 0) + (cumtrapz(rho(x)/eps_val(2))*dx).*(x > 0)];
E = @(x,eps_val) [(cumsum(rho(x)/eps_val(1))*dx).*(x <= 0) + (cumsum(rho(x)/eps_val(2))*dx).*(x > 0)];
eps_r = [10; 15];

figure(1)
plot(x, E(x,eps_r))
grid

The idea behind the ‘E’ function is that it does the logical discrimination (emulating the if block in your code) by multiplying by a logical (0,1) vector, essentially adding the two parts together, with one having a non-zero value for the appropriate range of ‘x’ for each part, and the other being zero for the same range. Together they produce a vector of non-zero values for the entire range of ‘x’. You can experiment with it by taking it apart to see how it works.

[Anton Velychko]

That works! Thank you very much indeed!

The only problem left now is that I have another couple of variables Phi and EE that I calculate from E using similar syntax to what you proposed for E, namely:

Vbi=0.27;
  Phi = @(x) -cumsum(E)*dx;
  EE1= @(x) (V_bi-Phi);%+EE_0;
  EE= @(x) EE1-EE1(length(x));

However, when I try to plot any of those:

   figure(1)
    plot(x, Phi(x))
    grid on

I get this error message:

Error using cumsum CUMSUM supports only numerical or logical inputs.

Could you advise how to rectify this?

Thank you

[Star Strider]

My pleasure!

I defined ‘E’ as a function, so cumsum will have indigestion if you give it the ‘E’ function as an argument. You could create ‘E’ as a vector by defining the vector separately:

Ev = E(x,eps_r);                        % Vector Representation Of ‘E’

and then cumsum (and the more accurate cumtrapz) will be happy.

If you have functions you want to integrate as functions without first evaluating them as vectors, I would use the integral function. It was introduced in R2012a, and there are other versions such as quadgk, and in earlier releases, quad, and quadl that work similarly. I would experiment with these, instead of evaluating them to vectors first. They’re adaptive numerical integration functions, and will be much more accurate than summing a vector. It depends on what you want to do.

[Anton Velychko]

Fantastic! That works brilliantly!! Many thanks!!! And cumtrapz is indeed a bit better than cumsum, but both giving an error of about a few % depending on step dx (the smaller the better obviously).

My eventual goal is to be able to plot all 3 variables - Ev, Phi and EE vs. x. With cumsum or cumtrapz that is easy because they both give results as a vector who's elements correlate with vector x [i.e. EE(1) = x(1) ... EE(n)<=> x(n)]. I appreciate that integral and quad give much better accuracy, but they output a constant value (correct me if I am wrong), whereas I need the answer for each of the 3 variables as a vector correlated with vector x or as a function of x.

If you could let me know how to achieve this with integral or quad - I would be very grateful.

Thank you

[Star Strider]

My pleasure!

Since you have integral, I’ll give it a go with it because it’s more powerful than its predecessors. If you have a reference outlining what you want to do, that would help. If it’s a paper, please direct my attention to the relevant parts and equations.

[Anton Velychko]

Hello Star Strider,

sorry for delay, it took me a while to write all up. Attached is a pdf file with geometry sketch, formulas and matlab code.

SUMMARY: To understand the math, you mostly need page 1 with the geometry and page 3 with the main formulas (10-12). Page 2 includes derivation of E and Phi within different layers, from which the total E and Phi are calculated. Pages 4 and 5 is matlab code I use with the current problems listed at the bottom of it.

IN DETAIL:

PAGE 1. There are 3 layers: A, W and D. A has eps1, while W & D - eps2 (same for both). In the figure - top (grey) is charge distribution, middle (BLUE) is Ev(x) and bottom (RED) is Phi (x), all qualitative sketches. Please ignore the formulas on the RHS, the main equations used (Eqs.10-12) are on page 3.

PAGE 2 Derivation of E and Phi components within the 3 layers (A, W and D).

PAGE 3 The main equations used in generic form (Eqs.10-12). Please note that all integrals are parametric (functions of x).

PAGES 4 & 5 Matlab code. The current problems I face are listed at the bottom of the code.

If you could have a go on adopting integral and quad - that would be fantastic.

Thanks a lot

[Star Strider]

Thank you. That’s far from my areas of expertise, so I’ll do my best to understand it and see if I can get it in a form that would use integral, since it is superior to numerically integrating vectors. I can’t promise anything, though.

[Anton Velychko]

Certainly, I do not expect any promises, you've already done me a great favour. :-) Please let me know if anything is unclear or confusing. Thank you.