3
$\begingroup$

I am trying to calculate an impedance value from voltage and current functions. The calculation method is based on how an actual measurement device performs it.

I have two different current function, where one is slightly more complex.

iBV[t_, freq_] = 
  i0*(-Exp[(alpha - 1)*bvVor*(U[t, freq] - E0)] + 
  Exp[alpha*bvVor*(U[t, freq] - E0)]) + Cdl*DU[t, freq];
i[t_, freq_] = iBV[t, freq]*iLimit/(Abs[iBV[t, freq]] + iLimit);

If i use iBV directly, the calculations finish very rapidly. When I however use i I am still waiting for the results after half an hour and I don't understand why. Granted, it has to calculate iBV twice for each round but shouldn't that lead to twice the time and not 1000 times?

As you'll see I need to integrate. Is maybe the complexity of performing the integration much more complicated in the latter case for some reason?

The rest of the computation is as follows:

ClearAll[t, freq, U0, Harmonic, bvVor, i0, alpha, E0, T, z, cBulk, 
  Cdl, iLimit];
ClearAll[U, DU, iVoigt, i, iL, Rea, Ima, tmax, ReZ, ImZ, Impe];

(* Current parameters *)
i0 = 0.0001; alpha = 0.5; E0 = 0; T = 293; z = 1; cBulk = 25; Cdl = \
0.000001; iLimit = 0.0005;

(* Voltage parameters and Misc *)
U0 = 0.1;
Harmonic = 0;
bvVor = z*96485/8.314/T;

(* Voltage Functions *)
U[t_, freq_] = U0*Sin[2*\[Pi]*freq*t];
DU[t_, freq_] = 2*freq*\[Pi]*U0*Cos[2*freq*\[Pi]*t];
(*iVoigt[t_,freq_]=Cdl*DU[t,freq]+U[t,freq]/100;*)

(* Current Function *)
iBV[t_, freq_] = 
  i0*(-Exp[(alpha - 1)*bvVor*(U[t, freq] - E0)] + 
      Exp[alpha*bvVor*(U[t, freq] - E0)]) + Cdl*DU[t, freq];
i[t_, freq_] = iBV[t, freq]*iLimit/(Abs[iBV[t, freq]] + iLimit);

(* Functions to calculate the real and imaginary part of the current \
by integration *)
ReaF[t_, freq_] = i[t, freq]*Sin[(Harmonic + 1)*2*\[Pi]*freq*t];
ImaF[t_, freq_] = i[t, freq]*Cos[(Harmonic + 1)*2*\[Pi]*freq*t];

(* Integration bound, 1 period *)
tmax[freq_] = 2/freq;

(* Real and Imaginary part of the current *)
ReI[freq_] = 2/tmax[freq]*Integrate[ReaF[t, freq], {t, 0, tmax[freq]}];
ImI[freq_] = 2/tmax[freq]*Integrate[ImaF[t, freq], {t, 0, tmax[freq]}];

(* The Impedance *)
Impe[freq_] = U0/(ReI[freq] + I*ImI[freq]);

(* Calculation of results *)
Freqs = 10^Range[-2,5,0.1];
Results = Impe[Freqs];

(* Output *)
Res = Re[Results];
Ims = Im[Results];
marker = Graphics[{Red, Disk[]}];
ListPlot[Transpose[{Res, -Ims}], {Joined -> True, 
  PlotMarkers -> {marker, .03}, PlotStyle -> {Orange, Thick}}]
TraditionalForm[TableForm[Transpose[{Freqs, Res, Ims}]]]

If you want to test the difference, just rename i[t_, freq_] to something, and rename iBV[t_, freq_] to i[t_, freq_].

I also tried rewriting the whole stuff in Function/Module nomenclature, getting rid of doubled calculations in the process, but that just made everything slower (i guess due to overhead<->actual evaluation time bias).

Edit

I tried to do some benchmarks:

Timing[Do[iBV[0.1, 100], {10000}]]
Timing[Do[i[0.1, 110], {10000}]]

This does indeed report what I thought: {0.328, Null}, {0.657, Null}. So each function evaluation takes twice as long for the more complex function.

I also further narrowed it down: Everything is fast up to the integration.

I changed - as Edmund suggested in the comments - all definitions to SetDelayed, which did not improve the speed though.

$\endgroup$
  • 1
    $\begingroup$ All your function definitions are using Set (=). You should be using SetDelayed (:=). For example U[t_, freq_] := U0*Sin[2*\[Pi]*freq*t];. Notice the colon :=. $\endgroup$ – Edmund Sep 10 '15 at 9:36
  • $\begingroup$ Ok, makes sense, I changed it. Did not change the problem, though :-) $\endgroup$ – Jens Sep 10 '15 at 10:10
  • $\begingroup$ Another, most likely unrelated tip: do not use capital letters for defining your own symbols (a symbol can be anything, variables, constants, functions...). Capital letters are used by system functions and thus using them can result in unpredictable behavior from time to time. $\endgroup$ – Sascha Sep 10 '15 at 10:21
  • $\begingroup$ Good point as well! $\endgroup$ – Jens Sep 10 '15 at 10:44
  • $\begingroup$ I'm not sure what "calculations" are being compared, esp. since the update indicates the calculations' timings are as expected, but you might try NIntegrate instead of Integrate. Or check if the symbolic integral can be found with exact coefficients -- for example, i0 = 1/10000 instead of i0 = 0.0001. And don't use SetDelayed instead of Set, at least for the definitions involving Integrate in that case. $\endgroup$ – Michael E2 Sep 10 '15 at 10:58
4
$\begingroup$

The source of your problem is using Integrate with i[t, freq] in the expression.

It doesn't integrate.

Using NIntegrate (as suggested by Michael E2) the problem can be solved in a reasonable amount of time. Below are two changes.

ReI[freq_] := 
 2/tmax[freq]*NIntegrate[ReaF[t, freq], {t, 0, tmax[freq]}]
ImI[freq_] := 
 2/tmax[freq]*NIntegrate[ImaF[t, freq], {t, 0, tmax[freq]}]

I also modified the function definitions by replacing = (Set) with := (SetDelayed).

I think you also should either give Impe the Listable attribute or use Map as shown below.

Results = Impe[#] & /@ Freqs

{332.822 - 0.00389736 I, 332.822 - 0.00490649 I, 
 332.822 - 0.0061769 I, 332.786 - 0.00777461 I, 
 332.822 - 0.00978973 I, 332.822 - 0.0123245 I, 332.822 - 0.0155157 I,
  332.822 - 0.0195331 I, 332.822 - 0.0245907 I, 332.822 - 0.0309578 I,
  332.822 - 0.0389736 I, 332.822 - 0.0490649 I, 332.786 - 0.0617559 I,
  332.822 - 0.0777626 I, 332.822 - 0.0978973 I, 332.822 - 0.123245 I, 
 332.822 - 0.155157 I, 332.822 - 0.195331 I, 332.821 - 0.245907 I, 
 332.821 - 0.309578 I, 332.821 - 0.389736 I, 332.821 - 0.490648 I, 
 332.821 - 0.617688 I, 332.82 - 0.777622 I, 332.819 - 0.978965 I, 
 332.818 - 1.23244 I, 332.816 - 1.55153 I, 332.813 - 1.95324 I, 
 332.808 - 2.45894 I, 332.8 - 3.09553 I, 332.787 - 3.89685 I, 
 332.766 - 4.90547 I, 332.734 - 6.17487 I, 332.683 - 7.77221 I, 
 332.602 - 9.78165 I, 332.474 - 12.3084 I, 332.271 - 15.4836 I, 
 331.95 - 19.4692 I, 331.443 - 24.4636 I, 330.646 - 30.7055 I, 
 329.394 - 38.474 I, 327.441 - 48.0797 I, 324.42 - 59.8389 I, 
 319.808 - 74.0175 I, 312.901 - 90.7301 I, 302.841 - 109.791 I, 
 288.734 - 130.543 I, 269.913 - 151.741 I, 246.3 - 171.642 I, 
 218.689 - 188.366 I, 188.743 - 200.469 I, 158.575 - 207.39 I, 
 130.158 - 209.537 I, 104.868 - 207.972 I, 83.3488 - 203.979 I, 
 65.645 - 198.712 I, 51.4232 - 193.038 I, 40.1789 - 187.519 I, 
 31.3747 - 182.47 I, 24.5166 - 178.029 I, 19.1847 - 174.226 I, 
 15.0388 - 171.028 I, 11.8104 - 168.374 I, 9.29146 - 166.193 I, 
 7.32143 - 164.412 I, 5.7772 - 162.966 I, 4.56416 - 161.797 I, 
 3.60948 - 160.855 I, 2.8569 - 160.097 I, 2.26279 - 159.489 I, 
 1.79325 - 159.001 I}

Then you can plot it

Res = Re[Results];
Ims = Im[Results];
marker = Graphics[{Red, Disk[]}];

ListPlot[Transpose[{Res, -Ims}], {Joined -> True, 
  PlotMarkers -> {marker, .03}, PlotStyle -> {Orange, Thick}}]

Mathematica graphics

and make the table

TraditionalForm[TableForm[Transpose[{Freqs, Res, Ims}]]]

Mathematica graphics

$\endgroup$
  • $\begingroup$ That works very well. I tried it as well based on Michaels comment and it was indeed the bottleneck of the calculation. Very nice additions as well, I will check the attributes out as I don't know them yet. $\endgroup$ – Jens Sep 10 '15 at 15:12
  • $\begingroup$ @Jack LaVigne I pasted Jens code and followed your suggestions, NIntegrate and Impe[#]&/@Freqs. I get the error message: NIntegrate::nlim: t = 2./freq is not a valid limit of integration and after 6 minutes your result. The clue is not only NIntegrate, but mapping Impe. My question: Have you this error message too? $\endgroup$ – user31001 Sep 10 '15 at 16:20
  • $\begingroup$ @Jack LaVigne to avoid this message I would recommend ReI[freq_?NumericQ] := and ImI[freq_?NumericQ] := $\endgroup$ – user31001 Sep 10 '15 at 17:06
  • $\begingroup$ @Willinski Sorry, my mistake. I also copied Jens code but changed all the function definitions by replacing = (Set) with := (SetDelayed). I neglected to indicate that in the answer. If you do that I think things should go smoothly. $\endgroup$ – Jack LaVigne Sep 11 '15 at 17:31
4
$\begingroup$

You can increase the performance of i considerably simply by making sure that IBV is evaluated only once.

i0 = 0.0001;
alpha = 0.5; 
E0 = 0; 
T = 293; 
z = 1; 
Cdl = 0.000001; 
iLimit = 0.0005;
U0 = 0.1;

bvVor = z*96485/8.314/T;

U[t_, freq_] := U0*Sin[2*π*freq*t]
DU[t_, freq_] := 2*freq*π*U0*Cos[2*freq*π*t]

iBV[t_, freq_] := 
  i0*(-Exp[(alpha - 1)*bvVor*(U[t, freq] - E0)] + 
     Exp[alpha*bvVor*(U[t, freq] - E0)]) + Cdl*DU[t, freq]

i[t_, freq_] :=
  Block[{ibv = iBV[t, freq]}, ibv*iLimit/(Abs[ibv] + iLimit)]

Timing[Do[iBV[0.1, 100], {10000}]][[1]]
0.207059
Timing[Do[i[0.1, 100], {10000}]][[1]]
0.239106

The same technique can be applied to improve overall performance.

iBV[t_, freq_] :=
  Block[{u, v},
    u = bvVor*U[t, freq] - E0;
    v = Exp[u alpha] (1 - Exp[-u]);
    i0 v + Cdl DU[t, freq]]

With this version of iBV, I get

Timing[Do[iBV[0.1, 100], {10000}]][[1]]
0.173164

and, of course, the performance of i improves, too.

Timing[Do[i[0.1, 100], {10000}]][[1]]
0.216504
$\endgroup$
  • $\begingroup$ Thank you for your input. I did not know about the Block construct. I actually had the idea during my tests and implemented it using Module constructs, but I discarded it at some point. I found that the bottleneck was indeed in the Integrate step, as described just now in the answer by Jack LaVigne. Your answer is however very useful as well to further tweak the performance! $\endgroup$ – Jens Sep 10 '15 at 15:10
  • $\begingroup$ @m_goldgerg Nice improvement to the timing. Curious as to why you selected Block rather than With? $\endgroup$ – Jack LaVigne Sep 10 '15 at 15:15
  • $\begingroup$ @JackLaVigne. Block is measurably faster than With. $\endgroup$ – m_goldberg Sep 10 '15 at 17:28
  • $\begingroup$ @m_golderg Thanks for the tip. I was blissfully unaware. $\endgroup$ – Jack LaVigne Sep 10 '15 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.