3
$\begingroup$

I have some rather complicated function (related to the capacitance of coplanar waveguides for specific geometric parameters, for those who are interested) which I've trimmed down to the following for simplicity:

SampleFun[w_, s_, g_, h2_] := Module[{a, b, c0, k2, k2p, k4, k4p, q2},
  a = s/2;
  b = s/2 + w;
  c0 = g + w + s/2;
  k2 = a/b*Sqrt[(1 - b^2/c0^2)/(1 - a^2/c0^2)];
  k2p = Sqrt[1 - k2^2];
  k4 = Sinh[Pi*a/2/h2]/Sinh[Pi*b/2/h2]*
    Sqrt[(1 - Sinh[Pi*b/2/h2]^2/Sinh[Pi*c0/2/h2]^2)/(
     1 - Sinh[Pi*a/2/h2]^2/Sinh[Pi*c0/2/h2]^2)];
  k4p = Sqrt[1 - k4^2];
  q2 = EllipticK[k4]/EllipticK[k4p]*1/2*EllipticK[k2p]/EllipticK[k2];
  Return[q2]]

The source of the equation is a combination of theory found in Gevorgian et al (1995) (fairly legible) and Garg et al (1979) around p. 376. In its entire form it will give the capacitance per unit length of a coplanar waveguide resonator with finite ground planes and a double layered substrate of finite thicknesses.

Now, since this function depends sensitively on the geometric input parameters, lets say that we can always take w = 1500*10^-9, s = 30*10^-6, g = 50*10^-6, and h2 is a number that goes from 0 (lets take 0.1*10^-9 as we can't divide by 0) to 400*10^-9 for the range that I am interested in.

If we look at the function in this range, then it behaves in a peculiar way; for low values there seems to be an error in the calculation, and it gives exactly 0.

bad plot

Does anyone understand what is going on here, and how I could go about fixing it? I believe the problem stems from the EllipticK[k4p] part, as the rest seems to evaluate nicely, but it might be a combination of factors so I didn't leave out more. It is quite problematic for my calculation, as h2 is approximately 100*10^-9 in the geometry I am considering, which falls exactly in this gap.

Note, this is for Mathematica 10.0.2, which now that I write it I should definitely update as it might be something that is fixed in future versions.

$\endgroup$
  • 2
    $\begingroup$ Look at similar issues, e.g. Optimization of ODE with respect to the initial condition. Elliptic functions may produce some imaginary perturbations for certain ranges of arguments where they should be real-valued. Then one should use e.g. Re. Nonetheless I find you've defined your function in a suboptimal way. $\endgroup$ – Artes Apr 12 '18 at 18:58
  • 1
    $\begingroup$ Simply use higher workingprecision Plot[SampleFun[1500 10^-9, 30 10^-6, 100 10^-3, x 10^-9], {x, 1, 500}, WorkingPrecision -> 25] $\endgroup$ – Akku14 Apr 12 '18 at 19:09
  • $\begingroup$ Using Re brings me to the internal precision limit however, and using Block[{$MaxExtraPrecision = 1000},..] doesn't solve that either. Moreover, WorkingPrecision in the plot is giving me that the precision of the argument function is less than 25. $\endgroup$ – user129412 Apr 12 '18 at 19:55
  • $\begingroup$ Looking at it more closely, I believe that k4p is very close to 1, which is what the function has trouble with. So I should probably repose the question as being how to evaluate the EllipticK very close to 1. I could do a series expansion, but the value won't always be close to 1.. $\endgroup$ – user129412 Apr 12 '18 at 20:44
  • $\begingroup$ @user129412 When evaluating your function the system didn't work correctly so I guess it might be improved, I haven't my computer at the moment so I cannot provide a better function. $\endgroup$ – Artes Apr 12 '18 at 21:18
7
$\begingroup$

It is unfortunate that the OP neglected to give a source for the formula (always give references if providing an unfamiliar formula!) in the question; I will nevertheless proceed with my guess, based on a cursory look at the provided code.

The heart of the issue seems to be about the OP neglecting to check if his/her preferred convention for elliptic integrals matches with Mathematica's preferred convention. I have ranted about this in this math.SE post, but suffice it to say that the OP's formula seems to be assuming the modulus $k$ to be the argument, while EllipticK[] needs the parameter $m$ as an argument.

Thus:

SampleFun[w_, s_, g_, h2_] := Module[{a, b, c0, m2, m2p, m4, m4p},
      a = s/2; b = s/2 + w; c0 = g + w + s/2;
      m2 = (a/b)^2 (1 - b^2/c0^2)/(1 - a^2/c0^2);
      m2p = 1 - m2;
      m4 = (Sinh[π a/2/h2]/Sinh[π b/2/h2])^2
           (1 - Sinh[π b/2/h2]^2/Sinh[π c0/2/h2]^2)/
           (1 - Sinh[π a/2/h2]^2/Sinh[π c0/2/h2]^2);
      m4p = 1 - m4;
      EllipticK[m4]/EllipticK[m4p]*1/2*EllipticK[m2p]/EllipticK[m2]]

Plot[SampleFun[15*^-7, 3*^-5, 1/10, x 1*^-9], {x, 1, 500}, WorkingPrecision -> 30]

plot


Here is the long-overdue follow-through. The key insight here is that one can use the relationship between the complete elliptic integral of the first kind and the arithmetic-geometric mean, after which one can then use the homogeneity relation for the AGM to stave off cancellations and bad scaling. Applying this to the specific problem here, we have the following stabilized routine:

SampleFun[w_, s_, g_, h2_] := Module[{a, b, c0, ab, ac0, bc0, sab, sac0, sbc0},
      a = s/2; b = s/2 + w; c0 = g + w + s/2;
      ab = a/b; ac0 = 1 - (a/c0)^2; bc0 = 1 - (b/c0)^2;
      sab = Sinh[π a/2/h2]/Sinh[π b/2/h2];
      sac0 = 1 - (Sinh[π a/2/h2]/Sinh[π c0/2/h2])^2;
      sbc0 = 1 - (Sinh[π b/2/h2]/Sinh[π c0/2/h2])^2;

      1/2 * ArithmeticGeometricMean[Sqrt[sac0], sab Sqrt[sbc0]]/
      ArithmeticGeometricMean[Sqrt[sac0], Sqrt[sac0 - sab^2 sbc0]] *
      ArithmeticGeometricMean[Sqrt[ac0], Sqrt[ac0 - ab^2 bc0]]/
      ArithmeticGeometricMean[Sqrt[ac0], ab Sqrt[bc0]]]

which is now much less prone to subtractive cancellation or overflow than the previous version. (Try using it in the Plot[] command above without the WorkingPrecision option set.)

| improve this answer | |
$\endgroup$
  • $\begingroup$ Based on your math.SE post and and the literature, I indeed agree that I made a mistake in using the parameter m instead of the modulus k. Even with that fix though, one still requires WorkingPrecision -> 30 to get the plotting to work. What would be the best way to have the formula give output by itself? Use something like N[m4p, 30]? $\endgroup$ – user129412 Apr 13 '18 at 11:52
  • $\begingroup$ It's backwards; Mathematica needs the parameter, while your references use the modulus. Yes, the WorkingPrecision -> 30 is a bit of a crutch; I have an idea on how to stabilize the formula, but you'll have to wait a while, since I'm in the middle of moving to a new place. $\endgroup$ – J. M.'s technical difficulties Apr 13 '18 at 11:54
  • $\begingroup$ That is perfectly fine for me, I'd be very interested once youre settled and have some time. For now, this will do though, I can finish my design! $\endgroup$ – user129412 Apr 13 '18 at 13:43
  • $\begingroup$ I'd still be interested, if you could find the time at some point. Hope moving was smooth! $\endgroup$ – user129412 Aug 7 '18 at 13:20
  • $\begingroup$ Apologies for the much delayed response; please have a look now. $\endgroup$ – J. M.'s technical difficulties Sep 24 '18 at 14:13

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.